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1. Introduction 

l-l- Physical Problem 

Helicopter rotors operating in high-speed flight 
encounter a number of important aerodynamic phe- 
nomena. Two major features that dominate the flow 
on the advancing side of the rotor disc exist. The 
first, key feature is the presence of t ransonic flow con- 
ditions. Transonic flow imposes major limitations on 
the high-speed performance of the rotor. These lim- 
itations manifest themselves in high vibration levels, 
power divergence, noise, and component fatigue. The 
second key feature is the presence of a large vortex 
system near the rotor. The vortex system is com- 
posed of a series of helical vortex filaments generated 
at the tip of the blade. The following blade, which 
may be experiencing transonic flow, frequently inter- 
acts with these vortices. 

The interaction between a rotor blade and the 
vortices from a preceding blade can have a large im- 
pact on the blade aerodynamic environment These 
blade-vortex interactions (BVIs) cause large changes 
in local pressure which can occur over short periods. 
The pressure changes exist within a flow field which 
is, in general, transonic, unsteady, viscous, and three- 
dimensional. The vortex passage, therefore, acts to 
modify an already complicated flow field. 

A rotor interacts with a vortex under a wide 
range of relative orientations. However, the essential 
physics can be illustrated by considering a rectangu- 
lar blade of infinite aspect ratio interacting with an 
infinite line vortex at an angle ft. Johnson (ref. 1) has 
shown that this problem is steady in a coordinate sys- 
tem whose origin travels with the intersection of the 
blade centerline and the projection of the free vortex 

on the blade. (See fig. 1.) The steady coordinate 
system is 

y — y — Afx cot, ft j ^ ^ 

The speed at which the origin travels is a function 
of the angle ft. When ft = tt/2, the vortex is per- 
pendicular to the blade and the speed of the inter- 
action point is zero. (See fig. 2.) For increasing 
values of ft, the speed of the interaction point in- 
creases, but the problem remains steady. For ft = tt 
there is no spanwise flow dependence, and the prob- 
lem is now two-dimensional. (See fig. 3). How- 
ever, the cost, of this two-dimensional simplification 
is that this problem is now intrinsically unsteady be- 
cause the speed of the interaction point, is infinite. 
The blade-vortex interaction may then be classified 
by the two limiting conditions defined by ft - t/ 2 
and ft = tt. The first condition (ft = ir/2) may be 




Axes moving along 
x span with vortex 

y'-y~ cot 0 



Figure 1. Interaction of infinite aspect ratio Made with i„fi- 
nite line vortex. 
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Figure 2. Low- S |>eed interne lion between rolor and vortex. 




0 - x; Mr = oo 

Figure 3. High-speed interaction between rotor and vortex. 

called a low-speed interaction (LSI) because this is a 
steady problem even in the original coordinate sys- 
tem. The second condition (ft = w) may be called 
a high-speed interaction (HSI) because this is an 
unsteady problem even in the transformed coordinate 



system. Both LSI and HS1 represent real interactions 
which can have significant effects on the rotor aero- 
dynamics. LSI, for instance, is the principal type o 
interaction which occurs during hovering flight. LS 
affects rotor power and low harmonic loading. Hb 
occurs during high-speed flight and descents and af- 
fects noise, vibrations, and the higher harmonics o 
loading. Furthermore, HSI contains all the^physics 

of LSI- therefore, the capability to solve for HSI con- 
tains the ability to solve LSI. The solution of HSI is 
the main study of this report. 

1.2. Physical Model 

The solution of HSI requires the computation of 
the time-varying surface pressures during the vor- 
tex passage. Because the angle between the vor- 
tex axis and the blade is zero, no spanwise flows 
are induced and the flow can be assumed to be two- 
dimensional. The distance lie tween the airfoil and 
the vortex is assumed to be large enough to assure 
a basically inviscid flow; that is, the vortex does not 
distort the airfoil boundary layer. Any shock waves 
present are assumed to be weak and not a source ot 
rotational flow. With these assumpt ions, the aerody- 
namic problem can be modeled by assuming a two- 
dimensional potential flow field. The mass conserva- 
tion equat ion for such a flow field is 

S+|.ij4 1 )+4(W=» <‘- 2 > 

at (>j- oy 

Since the flow is isentropic, the fluid density can be 
determined from t he Bernoulli equation 


P = 




l/(-/-t> 


(1.3) 


These two equations, solved together, constitute a 
full-potential model of the flow field. A solution for 
these equations can be accomplished b> using a hi c 
difference algorithm. Such an algorithm originally 
developed by Sieger and Caradonna (ref. 2) is pre- 
sented in section 2. This scheme is then modified to 
include various vortex models. 

1.3. Vortex Models 


A primary aim of this report is to study the effect 
of various vortex models on the solution of HSI within 
the framework of a 2-D potential finite-difference 
algorithm. All these models are candidates for use 


in 3-D methods. Four different methods have been 
used to model the vortex: 

1. The first model approximates the vortex effect 
as a change in airfoil angle of attack. The 
velocity field of a point vortex law is used 
to compute an induced velocity at the airfoil 
quarter chord. This induced velocity is then 
applied over the entire airfoil. This is referred 
to as the “angle-of-attack method." 

2. The second model is related to the first but, 
instead of imposing only a constant velocity on 
the airfoil, a distributed velocity field from the 
vortex is imposed on the airfoil surface. This 
is analogous to a lifting-surface method. 

3. The third model is to specify a branch-cut 
discontinuity in the potential field. The vortex 
is modeled as a jump in potential across the 
branch cut, the edge of which represents the 
center of the vortex. This is referred to as the 
“branch-cut method." 

4 The fourth method models the vortex by ex- 
pressing the potential as the sum of a known 
potential due to the vortex and an unknown 
potential due the airfoil. This is referred to as 
the “split-potential method.” 

The first two vortex models are typical of the linear 
integral flow methods which are used in all the cur- 
rently available rotor-analysis methods. Methods 3 
and 4 are flow models of the vortex and can only be 
used in finite-difference methods. An important as- 
pect of this work is to determine whether this more 
elaborate modeling is necessary. 

1.4. Historical Background 

The problem of blade-vortex interactions is cen- 
tral to helicopter aerodynamics because t he interac- 
tion of the rotor and its vortex system can have a 
large effect on the aerodynamic environment ol the 
blade. The four vortex models discussed in section 
1 3 reflect, the level of sophistication of the global 
theories within which they were developed. To ap- 
preciate fully the various vortex models, reviewing 
the basics of the global rotor computations is usetul. 

Because of the geometric complexity of the vor- 
tex system, early analysts (e.g., Glauert m 1926 
(ref. 3)) treated the wake influence on the rotor by 
using momentum theory and blade element approxi- 
mations. The resulting models led to simple algebraic 
equations for the induced velocities. Computing an 
effective angle of attack on each blade segment with 
the induced velocity at the disc is possible and is the 
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essence of the angle-of-attack method. The result- 
ing angle of attack is then used in conjunction with 
tabulated lift, drag, and moment data to compute 
rotor performance. This approach has proven to be 
useful in computing total aircraft performance but 
is unable to predict accurately the details of the ro- 
tor aerodynamic environment. Accurate prediction 
of these details could not be performed until the ad- 
vent, of high-speed computers (circa 1960). 

Piziali and Du Walt (ref. 4) proposed the first 
practical method for studying the details of the rotor- 
wake interaction. The complete rotor wake was 
modeled with a series of such line segments. The 
velocities induced by each segment at the rotor were 
combined to produce an effective angle of attack. The 
solution then proceeds as before. This improved the 
earlier model by allowing for individual blade-vortex 
encounters to be studied. Isay (ref. 5) presented a 

more general solution of the induced velocity for a 
spiral wake model. 

Kocurek (ref. 6) presented an extension of this 
method (for hover only) in which the blade was 
treated as a lifting surface. The induced velocity was 
computed over the entire surface of the blade not 
just at one point. These velocities were then used to 
compute local pressure and lift. Tabulated data were 
used to provide the drag and pitching moment with 
the solution proceeding as before. 

These approaches can be broadly classified as 
blade-element integral methods. Currently they are 
the most popular methods of rotor analysis in use 
(especially the angle-of-attack method). 'Numerous 
investigators have improved on the basic model cul- 
minating in the effort by Johnson (ref. 7). Despite an 
impressive versatility, however, the blade-element in- 
tegral methods have some serious shortcomings. Pri- 
marily the deficiencies of these methods lie in the 
use of tabulated airfoil data to provide aerodynamic 
forces and in the related assumption of local two- 
dimensional flow. Furthermore, unsteady aerody- 
namic effects are modeled with quasi-steady approx- 
imations which are incapable of modeling the truly 

unsteady phenomena of transonic flow. 

During the time period I960 1980, the field 
of computational fluid dynamics underwent rapid 
growth. Computer speeds increased to the point 
where it became possible to use finite-difference 
methods to compute simple rotor flows. Various in- 
vestigators (e.g., refs. 8 and 9) addressed the lim- 
itations of the integral methods by using finite- 
difference methods to compute rotor aerodynamics, 
typically with hybrid methods. The hybrid methods 


use finite-difference techniques to solve a limited part 
of the flow field and a linear int egral method to pro- 
vide the global .solution. This method provides the 
capability to compute the entire transonic nonlinear 
flow field near the rotor. An essential difference be- 
tween the nonlinear and the linear integral parts of 
the solution is that the linear solutions depend only 
on the blade surface and shear layer conditions be- 
cause the speed of sound is assumed to be constant. 
In contrast the nonlinear solution depends on flow- 
conditions in the entire flow field. For rotors, this 
field dependence is especially important because the 
held is frequently occupied with vortices from pre- 
vious blades. Therefore, an important part of the 
development of rotor finite-difference schemes is the 
means of specifying vortices. 

Caradonna, Tung, and Desopper (ref. 8) devel- 
oped the first finite-difference scheme that included 
vortices in the flow field. They solved a high-ti|e 
speed hover problem in which the vortices were spec- 
ified as edges of potential discontinuities (branch-cut 
method). This scheme produced good comparisons 
with pressure data. Interestingly, reference 8 also re- 
ported an inability to obtain a good solution when 
the effect of the vortex was included only by a blade 
surface inflow specification (angle-of-attack method) 
Strawn and Caradonna (ref. 9) solved a similar prole 
lem bv using a split-potential model. Their method is 
a modified version of a full-potential algorithm devel- 
oped by Bridgeinan, St.eger, and Caradonna (ref. 19). 

To date, the forward flight computations have re- 
lied on vortex-induced surface inflow boundary con- 
ditions (angle-of-attack method) and have been' fairly 
successful at high advance ratios where the induced 
flow is a small percentage of the total inflow. Never- 
theless, there remains a serious question of how best 
to introduce moving vortices into a computational 
grid and thereby predict their effect. The solution of 
the 2-D HSI is a convenient testing ground for the 
vortex modeling schemes that are required for the 
full 3-D problem. 

A number of investigators (refs. 11 to 17) have 
studied the 2-D HSI problem by using finite- 
difference methods. These investigators have been 
primarily interested in acoustic effects and have used 
modified versions of earlier algorithms and vortex 
models. The problem of vortex specification has not 
been a primary aim of these studies. George and 
Chang (ref. 11) modeled the vortex with the angle- 
of-attack method to investigate the effects of blade- 
vortex interactions on noise. Later they extended 
their methodology and results to reflect the results 
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from the branch-cut, method (ref. 12). McCroskey 
and Goorjian (ref. 13) introduced the split-potential 
method, which was first proposed by Steinhoff. (See 
ref. 14.) Computations with Euler and thin-layer 
Navier-Stokes algorithms have been presented in ref- 
erences 15 and 16. Sankar and Malone (ref. 17) pre- 
sented a full-potential solution by using a strongly 
implicit procedure. All these methods produce re- 
sults which are characteristically similar and, to 
a limited extent, show good correlation with each 
other, especially those methods that employ the split- 
potential vortex model. 

1.5. Purpose of Current Research 

Of all these methods, the full-potential approach 
is probably the best suited to rotor computations 
because it is geometrically general and is much faster 
than the Euler and Navier-Stokes methods, lhe 
purpose of the current research is to explore fully 
the BVI phenomena within the context, ol the tull- 
potential algorithm. This exploration requires four 

major tasks. 

The first task is to develop a full -potential algo- 
rithm which is specifically tailored to solve the BVI 
problem. Improvement to previous methods include 
(1) special boundary conditions to increase flexibil- 
ity in modeling airfoils, (2) a "full-potential" mesh, 
(3) allowance for a variable time step in the unsteady 
solution, and (4) an improved method of comput- 
ing the flux metric. Details of t hese improvements 
are presented in section 2, which includes a complete 
derivation of the full-potential algorithm. 

The second task is to implement the various meth- 
ods of modeling a vortex. This task can best be 
accomplished within the framework of a single algo- 
rithm. Using a single algorithm eliminates any ques- 
tion of differences raised by issues such as grid or 
boundary conditions. The present method provides 
a unique opportunity to accomplish this. The full- 
potential algorithm is modified to include the effects 
of four different ways to compute the influence of the 
vortex. These modifications include 

1. Modifications to airfoil boundary conditions 
for the angle-of-attack and lifting-surface 

methods 

2. Special internal boundary conditions to imple- 
ment the branch-cut method 

3. The inclusion of a split-potential model 

The split- potential method has not yet been fully 

implemented for unsteady BVI problems with the 


full-potential flow model; this will be accomplished 
for the first time in the present work. Details of the 
vortex modeling are presented in section 3. 

The third tasks is a unique side-by-side com- 
parison between the models. Comparisons with 
experimental data for subcritical and critical flow 
conditions are also made. These comparisons high- 
light the merits of the various models. This work 
provides the basis for a more systematic approach 
to three-dimensional computations of blade-vortex 
interaction. The results of these comparisons are pre- 
sented in sec tion 4. 

The fourth task is to make a parametric study of 
the BVI problem, which will be the first complete 
study presented in the literature. The parameters 
studied fall into two categories: ( 1) flow-field param- 
eters and (2) vortex parameters. The results of these 
parametric studies are presented in section 4. 


2. Full-Potential Algorithm 

The aerodynamic problem of HSI will be modeled 
with a potential flow-field assumption. Under this 
assumption the basic equations of fluid dynamics 
(mass, momentum, energy, and equation of state) 
are reduced to a system of two equations with two 
unknowns: the mass conservation equation, 

A(p*,)+ J-(P<&J/) =0 (2.1) 

dt dx dy 

and the Bernoulli equation, 

/, = f i + (w* - - *;)] 1/1 ' " 

L (2.2) 

In these equations, all velocities are normalized by 
0x ; distances, by the airfoil chord; time, by c/a ? c ; 
and density, by its free-stream value. 

2.1. Conservative Formulation of 
Transformed Equation 


The system of equations (eqs. (2.1) and (2.2)) 
is transformed to a computational plane under the 
general transformation 


£ = £(*m/T) ' 

T) = T)(x,y,t) > 

T-t 


(2.3) 
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and the equations must be conservative in these 
eoordmates. Equation (2.1) is written in conservati ve 
(or divergence) form; that is, 

a/+ ^ = ° M 

where F, is the flux of the quantity / (p for the 
potential model). This generic form must he main- 
tained through the entire .solution process (includ- 
ing discretization) if mass conservation is maintained 
V iv, and (ref. 18) has derived a general conservation 
form lor such a generalized coordinate system which 
transforms equation (2.1) into 


(taw/™ Vi 6 f, aCObian ° f the transformation 
(eqs. (2.d ). finder this transformation, the Bernoulli 
equation becomes 


P- 1 + 


M'i - 2$ r - U$c - V# 


l/h-l) 


where U and V are the contravariant velocity vec- 
tors, with U being the velocity perpendicular to the 
// direction and V being the velocity perpendicular to 

the £ direction. In general these velocities are defined 
as 

U = + Aj$r + A 1 

V = // r + A 2 $c-K4 ;! $ v J (2 ' () 

where Aj, A 2 , and A :i are metric terms defined as 

A, = £ + ^ 

A 2 = V{- Vi) = £,.//, + , (2.8) 

A; t = Vt ) 'v?}= + 

2.2. Computational Grid 

For the present purposes, three characteristics 
are useful for the computational grid to have The 
hrst characteristic is orthogonality, which is useful 
because it reduces the complexity of an algorithm 
The second characteristic is that the grid lines con- 
form closely to the shape of the airfoil; this increases 
the accuracy of the solution. The third character- 
istic ,s that the grid lines should align with the free 
stream away from the surface in order to facilitate the 
branch-cut method of modeling the vortex. For these 
reasons, an orthogonal H-mesh was chosen to be used 


with the current method. The streamlines and po- 
tent, at lines which surround an airfoil in incompress- 
ible flow form such a grid. This type of grid may 
be computed by means of a complex mapping solu- 

n. The Joukowski airfoil transformation is used in 

he present method because it provides a convenient 

closed-form solution. The grid is generated with the 
following steps: 

1 Produce a satisfactory stretched Cartesian 
grid using any method 

2. lise the (£,//) coordinates along the front face 

of the grid to integrate § to the aft face of the 

grid (this solves for the streamlines around a 
circle) 

T Transform the circle solution by using the 
Joukowski transformation to produce the air- 
foil solution 

T Select an appropriate distribution of points 
along the airfoil “streamline" 

b. Interpolate *(£.//) to find the potential at 
each of these points 

b. Find the location of each of the “off air- 
foil streamlines which have matching poten- 
tial values 

7. Form the grid with the result ing set of points 

The details of the development are presented in 
appendix A. Since the grid is orthogonal, the metric 
erm A 2 is identically zero. Furthermore, since the 
grid is steady, the contravariant velocities become 


U = A,$ f 
V = A;,*,, 

2.3. Metric Derivatives 


rile transformation of the equations to the com- 
putational plane gives rise to metric terms (ie f 
//.„) as a result of the chain rule. I„ matrix form ’ the 
chain rule expansion is 


J ’r j)r 1 

!/C o 
a 7/ I /»/ 0 


The determinant of th 


<>j I r d T 

Oy = 

(>! . d„ 


e matrix is 


(2.10) 


( 2 . 11 ) 





( 2 . 12 ) 


(’■ranter's rule can be used to solve for 

yj - ytJ 1 '! 1 

a = K 

'tV'4 
d >= 0„ - 

These expressions are applied tof and IJ, respectivelj, 
to obtain v 




(sy ) = - 


D, 


(2.13) 


(?,w) ~ D„ 

3n = D„' 

The terms </„, *. *„■ «d *< are called the primitive 
metrics and are detertnined by the follow, ng simple 
central difference formulas (for convenience A£ 

A,/= 1): 


, v _ ifi+i 

1) 


2 

_ (W + 1 

- 1 ) 

\y0ij 

2 

, , - iiiii 

- J 'j- 1) 

( J 0 ) / j — 

2 

1 , , . - { hil 

- J/j-l) 

( ytj )i j 

2 


(2.14) 


which exactly conforms to the airfoil, this leads to 
$ = 0. One problem with employing this boundary 

condition is that a new mesh must be generated with 
every new airfoil or airfoil orient ation. An alternative 
to computing a new grid is to use a transp. ratio., 
rather than a “noflow” boundary condition. 11ns 
approach uses a fixed grid which conforms to son^ 
convenient profile (e.g., a Joukowsk. airfoil) to ► ap- 
proximate the desired profile. The flow must th - 
fore pass through the grid surface at an angle ■ €, which 
is the angular difference between the grid surface and 
the actual surface. Figure 4 illustrates this relation- 
ship The flow normal to the grid surface * N is 


$ v tail f 


(2.15) 


where 4>s is the flow velocity tangent to the grid 
surface This condition is merely a generalization 
of the usual small- disturbance boundary condition. 
The value of to be used in the actual algorithm 
remains to be found. Since the coordinate system 
is orthogonal, the only difference between t he A and 



At the airfoil boundaries, the derivatives in the ij 
direc tion are computed with the aid of a pseuc ogri 
line inside the airfoil contour. This grid line is 
determined with a simple linear extrapolation of the 
grid points off the airfoil. 

2.4. Boundary Conditions 

The four boundary conditions that are imposed 
on the flow are associated with ( 1) the airfoil l surface, 
(•2) the outer boundary of the grid, (3) the aft face of 
the grid, and (4) the Kutta condition. 

I. Airfoil Surface 

For inviscid flow, the surface boundary condition 
requires that the flow be tangent to the airfoil sur- 
face This requirement can be met by setting the 
contravariant velocity vector V to zero. For a mesh 



(b) Flow -through condition on coordinate which approxi 
mates body (generalisation or smalbdisturbanceboundary 
condition). 

Figure 4. Transpiration boundary condition at airloil surlaie . 
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V directions is 
therefore, 


a simple one-dimensional stretching; 




? /A' 


(2.16) 


face so that. /> — 1 is ensured. With the Bernoulli 
condition, the following expression for is derived 
(see appendix B): 4 


Outer Boundaries 

Along the outer boundaries, the flow is required to 
return to undisturbed conditions. Often, computa- 
tional outer boundaries are so close to t he airfoil sur- 
ace that this condition cannot, produce an accurate 
or stable result . For these close-grid boundaries, spe- 
cial nonreflective boundary conditions are imposed. 
However, the current grid boundaries are sufficiently 
far away (155 chords horizontally and 80 chords ver- 
tically) that the assumption of undisturbed flow is 
valid The outer boundary conditions are set with a 
Dirichlet condit ion 


* = M x x 


(2.17) 


S.4-H. Aft Face 

Along the aft face of the mesh, the flow is also 
required to be undisturbed. However, because the 
present method employs a number of branch cuts 
(lines of potential discontinuity which model vortic- 
ity), the potential cannot be easily specified at this 
boundary. Instead free-stream conditions are im- 
posed by modifying the outgoing flux along the aft 





which is used in the flux computation. 


(2.18) 


2-4-4- hull a Condition 


for the lifting conditions, allowance must be made 
or a jump in potential across a wakelike branch cut 
(Kutta condition). This cut extends from the air- 
foil trailing edge to the alt, face of the mesh (which 
precludes the use of Dirichlet boundary conditions 
along the alt face). The cut is aligned with the mean 
chord line of the airfoil. In an unsteady flow, the 
jump m potential F across the cut must be convected 
downstream. The following equati on, governing the 
convection of vorticity from the trailing edge, is de- 
rived by using the Bernoulli equation and continuity 
ol density across the cut (see ap|>endix B): 

r r + (U)r f = () (2.19) 

where (U) is the average of the velocities above and 
below the branch cut. Equation (2.19) is used to 
determine the value of I along the branch cut. 


2.5. Difference Equation 


step* Differ ^ USef '! 1 111 ? eVe, ° pmg thB d'ffemice equation. The first criterion is tin 

of equations Th? 7 7 ■ i01 ' ^ # ^ ^ CriterioH *“***<* the „ee 

ot equations. The second criterion is that the con.servation form be maintained Tin 
accurate computation of shock motion and strength. maintained. Th, 

Equation (2.5) may be written by using a backward difference in time as 


at any current, time 
for solving a system 
criterion ensures an 


A + 1 

7-j -It-) +Af 


,A'+1 


puy' v +> f f>v\ ■ v + | 

Ff + (fc) 


= 0 


This difference equation maintain, the conservation form, but the density „i|| remain, „ 

11,15 Ca " be — «•* '» *'«■ d ™it} in a Taylor aerie., expansion ' ' + ‘’ 


; 2 . 20 ) 

level ; 


ri A -f 1 


dp 


N 


+ -*' v i 




( 2 . 21 ) 


The operator jfe is derived from the Bernoulli equation and is ex 


expressed in the following conventional form: 



1 dr + Of + eh) 


( 2 . 22 ) 
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A complete derivation o f e q u,ti ? («) i» 


t e v' + '_ (jl\ k = At-vf - (4-Y [*c+u*(c+v*„c) ; ' 

W VjJ VJ ..7 

/ ,J \ A” 1 q V— 1 

_ (— ) [# rT + U# {r + V* 1 ,r]' 

\ / 


(2.23) 


The terms in brackets in equation (2.23) are seen to have a form similar to the time denvatives m the following 
non conservative full-potential equation: 

* rr + 2 U# fT + 2 V* r<r = A, (p 7 ' 1 - Ai*?)*« + A:i(p 7 - 1 ~ A 3 *{)#„„ - 2 UV»,* ) 

Indeed, the term 


C* = A /> 


V iY ,}A'-1 


1 [$ TT + 


/V-l 


(2.25) 


,i u:« i pmi affects both the mass conservation 


( £\ N f ( * - v + * _ * « ) + (* v - * *- 1 ) + u a; (*? +i - *?) + v A ^ ; ' +1 

\Ju) 1 


+ A/ 


A'+ 1 


J, 


A’ + l 


•u) +^,(— V 

J // 


A +11 


(2.26) 


■ r , ..it ran | )( . computed at the N time step with an error of only 

it r:i: " ::: ;:rriL,: F „ r —— •« - » 

* (£e) V+l . <£ (*)' W" -•"»*< a-Y 

After applying the delta <0,1. dividing through by -<#**>. and collecting term*. e q u»t,on ,2.2,1) become* 

|l + AC V A| **«{+ A/ v A:ii„»,,- A, : 'Ac' +1 h (37)'’ 4 « 

-****" (ia-Mgy*’]'***'-**' 

- ($) ***** [a (x:) w + a (£) '"vi 

+ (2A1) (A^) ,*■'■- «*-■ + *■'-)+ 

A,- v [a, * £ v -'« £ + As< ,* s '1 (2 28) 




i;Y — 1 ' 


V p . 

Equation (2.28) may be written in the following compact form: 

W *A>t -«•'■)= *■(/,#») + (#*' - **-' > + C-,p' *^* •»*’'> 


(2.29) 
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where 


»' = (jfi I A<W+‘ 


and 


h ( yM s ( 9* + s„ 

«T( 


r) «***'] 


(2.30) 


C* = 


. -f-'V — i 


ti 


At 


/V 




At v -l 


(• A ' -2* A ‘ -1 + (j» a ’ -2 


> + ( yr) A/A (r A - 2 1 * : 


+ 


,fV 


A# a (A,$ a ' + A :^, 


■ f; <’r/ 


(#A 


(2.31) 


The left_hand side term L is simply the time update term which arises from the delta form of the equation The 

2.6. Implementing the Algorithm 

2 . 6 . 1 . Basic Difference Operators 

Equation (2.28) must now be implemented. The basic difference operator for time is 


= 


$A' + 1 


fc-V 


= *r 


(2.32) 


A/ v - 1 

The velocity term which is used in the computation of L and R. is formed with a central difference 

(«|> c ). - ^' +l ~ 1 

’ * h ~ 2 (2.33) 

The computation of the corresponding #„ term must allow for jumps in <t> across branch cuts. For convenience 
the gnd is assumed to Ire permeated with horizontal branch cuts which lie slightly above the grid mints Farh 
grid point, has a jump in potential F associated with it. The value of r is L2 “ , g P ?' Eadl 
actual branch cut. A difference expression which accounts for this field of branch cuts is ^ 


" (*j + fj) + *, -(#;_, + T ,-_,)] 


(2.34) 


figure 5 illustrates the velocity computation near a vortex branch cut. 

2.6.2. Flux Operator Differencing 

Each of the spatial flux terms is made up of the product of three terms: (1) a velocity term (2) a density 
gnu nodes, the velocity term is computed by using a one-sided difference 

(**)i+(l/2) = *('+1“ */ (2.35) 

BecauseVhe^fSI iZhl^ -T (swUrbtn f- The density terms are used to aid the stability of the algorithm 
■ — problem is transonic, the type of the equation will change from elliptic to hyperbolic depending 

i : H aC number - S : ability is achipved by -itching the type of operator with the equaZi tT To 

. trate the requirement for switching, consider the following two-dimensional equation: * 


1 - Ml) + 


* 


yy 


(2.36) 
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H 2 



Branch cut 

• r — ©; 

1 v v 


Figure 5. Computation of velocity at grid point in field permeated with branch cute. 


We define operators such that 


= 




(2.37a) 


A, # , = 


*«•+! -** 


(2.37b) 


It. has long been recognized that t he scheme 

(l - M’i) V j Aj 4» 4 Vi/A w 4» ■ 0 


is stable for M^c < 1, and that. 


(l - Ml) V,V,*4- VyA„* = 


(2.39) 


is stable for M x > 1. A problem arises when we try to use either of these schemes alone for a transonic flow. 

if < 1 and equation ,2.38) is convergent fo, M« > 1 only tf 


(\-M‘i)Ay\- 

,r iif , i Qf-Kilitv is achieved bv switching from equation (2.38) to equation (2.39) 

^ jJZSd that in. evaluated at the upstream cell center. (See sketch A.) A parameter rs then en, ployed 
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Spatial differencing and stablility: 

Stable difference for spatial flux in subsonic flow is 

Mx> x = V (Pi+ 1 /2 C> * — x ^ 


i + 1 

— • 


" (p i+l/2"5(Pi + Pi + l) 


Centered difference: 


Stability in supersonic regions requires upwing biasing and 
a simple way to accomplish this is 


(P^ - V(p i+1/2 A<1>) j^> i j 


Pi-3/2 ^Pi-1/2 


This remains a centered scheme; p is shifted however 
Sketch A 


!r ^ - ** ■ ^ .»* ^ 




( ] _ . . . \ P[ + I P± . Pi T Pi— I , 

(1 C) - + Ip 1 (#, + 1 - #,-) 


2 ' ' 2 
\Pi + Pi- 1 


(1 - ,,,_,) AlZiz. i+ Pi^i + Pi-2 


(#;- $;_j ) 


(2.40) 


The switching parameter /' is defined as 

d<e<io) 

u = 0 (is* < 0 (subsonic)) 

is= 1 


(2.41) 


(is >1 (supersonic)) 

The parameter C is used to provide additional numerical viscosity in the supersonic region The form of C is 
completeK. arbttrar^ In the present method, (7 is varied linearly with the La. Mach\umber by 


the following equation: 


C = 16M/ — 10.8 


by " u " , " ic '" experi,mni ^ criticai d ™ si,> - '• - *'-"""-<1 «■" tv 


2 / 7-1 \1 */(7-l) 

7+ i V + T“ Ux 


2.6. 3.2. Flux metric term (consistent differencing). Computation of the flux 
problem. The flux metric term is 


Ati - (fi + tr 


n / /+ 1/2 


J>+ 1/2 


metric terms poses a special 
(2.42) 


11 



The simplest, method of evaluating this is to average the values at the nodes so that- 

(AA \JL1l (2.43) 

This evaluat ion can lead to an error for a grid with stretching. The error occurs because the values 6 and % 
are computed by using central differences (eqs. (2.13) and (2.14)), and therefore, information from node (i- 1) 
to (i + 2) is used to compute the metric term at (i + 1/2). This “extra” data acts to reduce the accuracy ol 
the metric computation especially for regions where the metrics change rapidly. The efTect, of the diminished 
accuracy is to introduce an artificial mass into the flux computation. The extra mass can be removed by the 
use of a free-stream subtraction mat rix. This matrix is generated by specifying free-stream conditions on the 
mesh (4? = A / x j ) and computing the value of R x = R(4* = M^x). The matrix R x is then subtracted from 
R in subsequent computations to restore the mass balance. 

P’lores et al. (ref. 20) recently proposed a superior method of computing the metric which eliminates this 
problem. Their method involves calculating the flux metric terms at the same points in space at which the 
flux is differenced. The method is referred to as the “consistent-metric method.” Consider, for example, the 
incompressible mass conservation equation which must hold for the free-stream subtraction to be zero 

/Ai$A + = o (2.44) 

V Jji / c V / 1 / 

Here, the density is set to 1. and from the chain rule, 

#£ = 7-($x - 1U*T)) 

S J- 

>h 


—{&! , - >hl®v) 

7 ly 


(2.45) 


With equations (2.45) and (2.13), equation (2.44) becomes, for free-stream conditions, 

-1 f X ( V Ijt y^i ;) h 


(2.46) 


It therefore follows that the metric difference operators must commute in order to produce a zero free-stream 
subtraction. For this method, the primitive metrics (in the A’ direction) are computed by using a backward 

difference 


the 1 1 derivatives arc 


_ 1_ 

'/.+ I/-2J “ 2 




**.•+1/2.; “ * i+1 Xi 


y^M-u y>+ 1 ^ 

j j+i - x - 


(2.47) 




/ Vj ± 1 ~ Vj - 1 


i + 1 


j' f/j + l - Vj - 1 


(2.48) 


The determinant D,„ is 


7 ?. n — iftj 


(2.49) 
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and 


6 = 


y,, 

D,„ 


ih = - 


(« = - 


if 

A, 


(A 

V A, 


? >.v = 


D, 


Tlic flux metric term becomes 


J in — A 


A, 

J 


(2.50) 


'e+% 


The metric term ($f) ;+]/2 is computed in an analogous fashi 


(2.51) 


ion. 


;:;r , :r or -r — - -« — <«. - 

and one at the nodes. Although t he m^i i ^ °“ *“*' f °' thp * 311 d * A- terms 

coding, it provides a more robust method of "" , 101 ,llore storage space and complicates the 

method. The final form of t Z Tux ten„t ” ^ ,net " C8 ^ ^ been porated in the present 




' ' /- 


, , , Pi + Pi- 1 1 / A | \ 

+ "'~ (-J , <*f+l -*/l 

J \'*n J i + lft 


(l - 


v- 1 + p t - > 1 / A| 

2 JUJ,:- 


1/2 


(*/-*/- 1 


(2.52) 


«ju£^ion L ( 2.2$)^ i2A1) ' (2A0) ' (2 ' 35) ’ (2:M) ’ (2 - 3;5) ’ and < 2 "' !2 ) «re u*d in the implementation of 

2.6.3. Approximate Factorization 

A "".compact matrix invention a., follows is still required to solve equation {2 .2S): 

LA# V+1 = R v 


becomS^ "" l ° CO,, " >8 " r ° n "' lh - °l ,etator L » approximately factored. Tlte e„„a, 


Lion 


L,.LtA#' v+1 = R v 


(2.53) 


I'f T^' r '» l «o -'hit. at, error 

solution, and (:i) the overall scheme is stable. Tl . p°lU 3t h fc‘l"‘“T7 ° blai " *• 

two space derivatives, and this lead, to an ADI type T^faX L T ^ «“ * 


= [l+ A«*A 3 *,% - Af*A,* + , (ie) ft, (±i) />,] 


(2.54) 
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'I'lie form of Lc is similar. The final form of equat ion (2.29) is 

|i + a.**,*,*, - (£) «. (ft) 

x [i + A,* A - A,' Ae- (2J.) e f (ft) /A] <*" +1 - 


= Af A At’ v+l 




(ft) 


+ ,#*-•*-)+ l^srl Is 


j 9 a ’-‘\ / A/ v \ ,*.v _ , *.v-2 


£_ j ( $ v - *2$ V 4- $ ‘ 


+ 




..V- i ' 




n A a 4 A' / «A r _ ^ A : - 1 ) 


A/' V [Ai*f- 1 «£+ 1)+ V/^V ^ (/? 9 

Both operators Lc and L„ yield tridiagonal matrices. For example, 

L,,$* = A$*_i + B#* + C$* + i 


( 2 .55) 


(2.56) 


where 


A = 


B = 


-A r a 




pj + Pj - 1 






+A/ a a/ 


ii A 7VJJi-i/2 2 


c = 


A : A A 


Ac A3$ 




(2.57) 


^.6.^. Implementing Boundary Conditions 
the airfoil, the operator L, ( becomes 


L,,** = ** + 


A/‘ V+ 1 A:{<^i/A5 m At A At A j 1 \ ( A*\ £z+I±£-(* , +1 _*,)-( ^i) ft?AS„ 

[\J,,/j+(l/2) 2 A " 


d v 


(2.58) 


H the term AS refers to the change in velocity on the upper surface from time N to N + 1. Because this 
• Th^lTq^y " brought* to the right-hand side. The equation then becomes 


, R v-i a .a a ,^- v - 1 — 2 _ ( 2A A/ a A/‘ v+1 (Ai) PsS» 

A#*_1 + B*J + C** + , = R+ Ai At Aj« f , p S \ j \ 3 n ) s 


-Al V A:,*„Asf - Ar' A)' + l (^)(ft),^' 


(2.59) 
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where 


A = 0 


B = 


1 + 


At n + l At* 

0N 


£l\ Pj+l ± Pi 


C = - 


At v_rl A/ v 

( A -A 

8 A ' 

UJ 


J+l/2 2 

Pj + l + Pj 


2 


“ r;: rr^. “t* s-r - : he ^ - «* . «*. 

and can be determined by symmetry. P SUr ;U < ls s,ni,lar to that on 1 l|p upper surface 

2.6.5, Solution Steps 

The solution is obtained in three steps by the following equations: 

R 

L,A$ ;V+1 - 


L,,& 


(2.60) 

(2.61) 

(2.62) 


# A ' +1 =#A'+A#A'+1 

time step. Given tlus’Jduu^ Ztiri'l^ T ^ airfoil ^ 06x1 

" w<w - ^ *» * — -w 1 n ' von ” 

2.7. Steady-State Algorithm 

«|i.atio„ (2.55) i, "Tu" ° f 1 S,,!a<lj <' K>h,nn For ««*hr problems, 

tern. s which , he s„l„, ion. E(lu i tio „ " |Ua ' io " *>"»«» l»™<b ,i„„ 


1 1 - a/ v ( h. 


—ii 

8 -v 


/T v 

A/ v 


On 


A;i ) .v 0 
~r p (> ’i 

n 


J 1! 




J H 


8 A 

7l) /•%**] 


°*3f] 


(* 


A’+l 


-$ ;V ) 


(2.63) 


IS reduced to make the most efficient use nf M , ' ? ft A/ ' the error for a wide range of frequencies 

«e required ^^g~52 T T ^ «» 

the full algorithm (with all the , im e ,e„ ns i„ c |„*,“ ’ neptMe whe - **"*•> »00 atop, are required with 

The residual is determined by the following steps: 

J - (: °mpute the value of R at the first step 

2. Survey R to obtain its maximum local value Rj 

3. At each subsequent computational step, obtain R„ (the maxim,,,,, local value of R) 

4 ■ Drt « m " < ‘ v “'"- of K../R, Which is the normalized residual 

Wfc “ ll ' iS ralll< * reacll ' s 10-1 ‘he computation has converged. 
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3. Vortex Models 

In section 2. an algorithm for solving for the 
potential flow field around an airfoil in transonic 
flow was presented. This algorithm is now modified 
to include the effect of a two-dimensional vortex 
passing near the airfoil. The four models for the 
vortex discussed in section 1.3 are used: the an^e- 
of- at tack me thod, the lifting-surface method, the 
branch-cut method, and the split-potential method. 
These methods may he grouped into two categories. 
(See fig. (i.) 

The first category called the surface-specification 
methods model the effect of the vortex as an imposed 
normal velocity distribution on the airfoil surfac«_ 
Both the an gle-of- attack and lifting-surface methods 
fall into this category. These methods ongmate 
within linear-integral rotor theories. I he effect of 
the vortex on the general flow field is usually no 
considered in these- theories. These models are valid 
for linear flow fields. 

For problems characterized by the transonic non- 
linearity (that is, with a speed of sound that varies 


throughout the flow field), a surface effect cannot 
completely model the effect of the vortex; therefore, 
it is necessary to insert the vortex explicitly into the 
grid. This category is called explicit modf s. bo i 
the branch-cut and split-potential models fall into 
t his category. 

The vortex modeling liegjns with the ideal two- 
dimensional vortex potential: 


27T 


(3.1) 


where 9 is the angle subtended by the vortex and the 
field point. The tangential velocity at the field point 

is r 

r _ J_ (3.2 

1(1 - 

The singularity at r = 0 is the source of numerical 
instabilities and requires the use of an artificial I core. 
In the following tasks, the model developed bv Scull} 
(ref 21) given in the following equation is used: 


V o — 


2tt r 


r't + a.^ 


(3.3) 





Vortex potential 
or velocity field 
specification 


Blade inflow 
specification 


Explicit models 


Surface-specification model 


Figure G. Principal vortex models. 
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where a is the vortex “core” radius. The word 
core is used in the conventional sense that is used in 
the helicopter industry; that is, it refers to a region 
of rotational flow near the vortex center of rotat ion. 
This rotational flow may be restricted to a discrete 
region in some models or it may be modeled by a 
decay function as in equation (3.3). The radius" a 
defines the region in which the flow is rotational: 
within this region, the potential equation model of 
the vortex is invalid. 

The vortex is moved through the computational 
grid by integrating the flow velocity at the vortex 
over the current time step 


A-; V+1 = A ;. V + U v At N ) 
>« A+l = >-; v + v v Ai* ) 


(3.4a) 


The vortex convection velocities V,. and V,. can be 
determined by three different methods: (1) a priori 
specification, (2) interpolation, anti (3) the velocity 
nelcl of a point vortex. 


By far, the easiest method is to specify an initial 
position and then allow the vortex to convert at the 

free-stream speed. The equations most often used 
are 



Equations (3.4b) produce a “fixed-path” interaction 
Although specifying the vortex is a trivial matter 
m a 2-D flow problem, it is the usual procedure 
m 3-D integral computations of advancing rotors 
because of the cost and complexity of finding the 
wake deformation. This method is also very useful 
for various comparison studies. 


3.1. The Surface-Specification Methods 

The surface-specification models are produced by 
modifying the airfoil surface boundary condition^ 
based on an assumed vortex velocity field. This 
approach is basically the same as assuming that the 
vortex velocity may be superimposed on the general 
flow problem in the same way as a gust velocity 
would be modeled. The accuracy of this assumption 
depends on the location of the vortex with respect 
to the airfoil. If the vortex is far enough away, the 
held which it produces does resemble a gust field A 
constant velocity field produces the angle-of-attack 
method, and a variable velocity field produces the 
lifting-surface method. Furthermore, if the vortex is 
far enough away from the airfoil, the signal arrives at 
the various points along the airfoil at approximately 


the same time; such a solution would be a close 
approximation to the exact, solution of the flow field. 
However, the effect of a finite signal propagation 
speed is still violated by these methods and this 
becomes increasingly important as the vortex nears 
the airfoil particularly for transonic flow conditions. 
I he surface specification models are not capable of 
modeling the ‘time lag” between a signal arriving at 
a point on the airfoil surface near the airfoil and the 
signal arriving on the opposite side of the airfoil. 

3.1.1. Angle -of- Attack Method 

The angle-of-attack method is the simples! pos- 
sible model of the effect of a vortex on an airfoil. 
Equation (3.3) is used to compute the velocity at. the 
airfoil quarter-chord. With this velocity, a "vortex- 
induced angle of attack is computed. (See fig 7 ) 
The velocity perpendicular to the chord line is (if 
the leading edge is at r =z 0) 


...n-I ( 0-25 -A -„ 
-f a/ 


(M) 


The vortex-induced angle of attack is 


This angle is added to the airfoil angle of attack. The 
potential field is then computed as before. The angle 
of attack is updated at each time step as the vortex 
moves by the airfoil. 

I he airfoil is therefore assumed to lie a point in 
space. In order for this solution to be valid, the vor- 
tex must be far enough away for this approximation 
to be accurate (e.g the signal must arrive at every 
point on the airfoil simultaneously). 



Figure 7. A ngle-ol-atlack vortex model. 
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3.1.2. Lifting-Surface Method 

The lifting-surface method is the most general 
form of the surface-specification models. As with 
the ail gle-of-at tack method, equation (3.3) is used to 
compute the induced velocity at the airfoil surface. 
(See fig. 8.) However, unlike the angle-of-attack 
method, the velocity is allowed to vary over the 
surface. This provides a more physically accurate 
model of the effect. For the lifting-surface method, 
equation (3.5) is modified and is 


[ r 2 + a' 2 J 


' I = l 0 COS0 = — 


With this change, the computation proceeds as in the 
angle-of-attack method. 

The airfoil is therefore assumed to be a lifting 
surface. In order for this solution to be accurate, 
the signal must arrive at the upper and lower surface 
simultaneously. 


3.2. Explicit Models 


3.2.1. Branch-Cut Method 

Caradonna (ref. 22) was the first to use an explicit 
method vortex model with a finite-difference rotor 
computation. The branch-cut method, which he used 
for steady 3-D flows, is based upon the known poten- 
tial solution for a two-dimensional vortex (eq. (3.1). 
This potent ial is implemented by means of a branch 
cut. which extends from the center of the vortex to 
the aft face of the computational grid. (See fig. 9.) 
A jump in potential equal to I' is imposed across the 
cut. Because equation (2.34) has already been im- 
plemented to account for the wake cut, no changes 
in the algorithm are required. 



Figuie 8. Lifting-surface vortex model. 




«V r v 

•*/ 


Figure 9. Branch-cut vortex model. 


At first the branch-cut method seems to be well 
suited to a potential finite-difference algorithm. Dif- 
ficulties arise, however, in unsteady problems when- 
ever the vortex is moved. As the edge of the cut 
moves past a node, an abrupt change in the local 
potential occurs. This sharp change causes spurious 
waves which affect the entire flow field. The problem 
can be solved by spreading the edge of the branch 
cut, that is, by distributing vorticity on the vari- 
ous nodes which surround the vortex center. The 
simplest distribution involves the use of the near- 
est, four grid points. The distribution is weighted 
so that, the “center of gravity” of the vorticity repre- 
sents the center of the vortex. With four grid points, 
this will uniquely determine the vorticity distribu- 
tion. Increasing the number of points w ould require 
an arbitrary distribution to be imposed upon the vor- 
ticity. With this modification, the vortex may be 
moved from cell to cell smoothly and the spurious 
waves are reduced (not eliminated). The method of 
distributing the vorticity is illustrated in figure 10. 
The vertical distribution of vorticity is 




(3.8) 


and 


r 2 = r r -ri 


The horizontal distribution of vorticity is 


(3.9) 


x 1 — x v 

r 4 = i A i — 

X * 


( 3 . 10 ) 


and 


r-i = r 


X'2 


- Xv 


(3.11) 


The main vortex at (x r ,y v ) is then modeled with two 
branch cuts of varying strength; this can be called a 
two-cut model. The vortex may, in fact, be modeled 
by any arbitrary distribution. Stremel (ref. 23) uses 
a method in which the vortex is modeled with an 
area-weighted distribution of vorticity. A parabolic 
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Figure 10. Distribution of vort icily to four points nearest 
vortex. 


distribution of vorticity in the horizontal direction 
coupled with a linear vertical variation was used in 
reference 24. The effect of distributing the branch 
cuts is to create an artificial core for the vortex. The 
efficacy of the core is dependent upon the distribution 
of the nodes which are in the vicinity of the vortex. 
Because each of the separate branch cuts represents 
a separate subvortex and each subvortex has its own 
singular point, the core is very sensitive to the grid 
geometry. 

Another problem associated with moving the vor- 
tex is concerned with computing the vortex convec- 
tion velocity. Interpolation of the local velocities near 
the vortex is the only available means of comput- 
ing the vortex velocity directly. The interpolation is 
complicated by the fact that the vortex creates such 
a large local disturbance; separating the effect of the 
flow field on the vortex from the effect of the vortex 
on the flow field is difficult. One major shortcoming 
of the branch-cut. method is no good wav exists to 
separate the effects because the branch cut contains 
the combined potential of the vortex, t he free stream, 
and the airfoil. The vortex -induced velocities domi- 
nate the flow near the vortex and make an accurate 
interpolation very difficult. The velocity at the vor- 
tex may be approximated by using the velocity field 
of a point, vortex in conjunction wit h the lift on the 
airfoil. In t his approach, equation (3.3) is essentially 
used in “reverse" with f being the jump in potential 


at the airfoil trailing edge (and hence a measure of 
the airfoil lift). 

Another feature of the branch-cut method is the 
fact that it, requires a difference equation (eq. ( 2 . 34 )) 
to implement the effect of the vortex. That is, the 
vortex effect, is specified entirely by the potential 
jump which is represented by t he differencing across 
the branch cut. The accuracy of this difference is also 
dependent upon the local grid geometry. Therefore 
the accuracy of the vortex model changes as the 
vortex moves through the mesh. The distribution of 
the vorticity on the mesh further distorts the model 
bv increasing the mesh dependence. A successful 
branch-cut model is therefore a compromise between 
an effective core model and an accurate vortex model . 

3.2.2. SpUl-Potential Method 

An alternative to the branch-cut method is the 
split-potential method. (See fig. 11.) I n this ap- 
proach, the velocity is assumed to be a combination 

of a known velocity and a. perturbation velocity as 
follows: 

q = V6 + V g (3.12) 

where V G is the known velocity field solut ion and Vd> 
is a perturbation velocity, which need not be small. 
The total potential for the blade-vortex problem can 
be split between the perturbation potential (associ- 
ated with the airfoil) and the potential G, which de- 
scribes the vortex velocity field, as follows: 

* = o+G (3.13) 

Any potential algorithm may be modified in this way 
to include the effects of a known velocity component 
and a. perturbation velocity. Furthermore, the po- 
tential G need not represent a. vortex but can in fact 
represent any flow field which independently satis- 
fies the potential equation. When equation (3.13) 



Figure 11. Splil-jxiteiitial vortex model. 
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is applied to equation (2.28), the following equation 
results (in the compact form similar toeq. (2.29)): 


L(/.Ao v+, ) = R*(C V .« V +G V ) 



( v ,d Y_l ) 

4- C’(G v , G v_l , G v_2 , ij N , v - 

') 

+ Ao v f AG V -L(/- V .AG V+I ) 

(3.14) 

The Bernoulli equation undergoes a 
modification 

similar 

= p{<t> S + G A ) 

(3.15) 


The left-hand side of equation (3.14) is identical 
to the original algorithm (eq. (2.28)). The right- 
hand side contains additional spatial and temporal 
gradient terms in G, including the update operator 
L(G). 

Implementation of equation (3.14) has proven 
to be a challenge to several researchers whojiave 
sought to simplify the equation. (See refs. 7, 13, 
and 17.) The principal focus of these studies has 
been in eliminating the temporal gradient terms in 
G. These terms pose a particular problem because 
they involve the potential G explicitly. ( 'Om pitting 
these terms requires the tracking of a branch cut 
through the flow, in effect, the split-potential model 
is reduced to a branch-cut. model. This will be 
particularly difficult for the complex geometry of the 
full 3-D problem. Furt hermore, it is advantageous to 
minimize the computational requirements as much 
as possible. McCroskey and Goorjian (ref. 13) have 
shown (for a small disturbance formulation) that 
L(G) and AG V can lie eliminated since the vortex 
potent ial (eq. (3.1)) is a solution to 

LAG = 0 (3.16) 

Therefore, a "small disturbance version” of equa- 
tion (3.14) would be 

L(/\A<3 A ' +1 ) = RV V , 4> S +G V ) + (3.17) 

because C* terms are not present in a small distur- 
bance form. Sankar and Malone (ref. 17) restricted 
the solution to his algorithm to a so-called ‘weak 
split-potential’' approach in which the temporal gra- 
dient and most of the spat ial terms in G are simply 
dropped. 

The present method is neither asmall disturbance 
or weak split-potential form. In spite of this, the 
algorithm can be simplified by using a method pro- 
posed by Hoger Strawn, Aeroflight dynamics Direc- 
torate. U S. Army Aviat ion and Missile Command, 


Ames Research Center. The first step is to recast 
equation (3.14) as 

L(e v . A* v+I ) = R‘(/- V .c- V + G v ) 

+ cV, /"'>*. .** 0*"') 

+ C*(G ;V , G v_l , G* - *,##*, d v_l ) 

4 Ao v + AG v (318) 

The update operators L have been recombined to 
include the total potential. The temporal conser- 
vation correction term C* in G has been retained, 
since there is no way to effectively separate the G 
and 0 parts of the density. The effect of these 
terms on the solution is discussed in section 4. So- 
lution of equation (3.18) is achieved as before with 
equations (2.60), (2.61), and (2.62), which produces 
<j> ; Y+ 1 \ final step is added to obtain <f> 


■ - AG‘ V+1 (3.19) 


In implementing these equations, the velocity due 
to the vortex is computed by equation (3.3). The 
velocity components are 


G x — 


L ( y-V" \ 

2i r \r 2 +a 2 / 


(3/20) 


and 



( 3 . 21 ) 


These equations are transformed into G^ and Gy by 


Gc = GjrXt + GyPt (3.22) 

and 

G t , = GxXi) + G y y,, (3.23) 

The time derivative terms are obtained with the 
chain rule as 


G r = G^ r + GjjJir + G// r (3.24) 


This value may be determined most easily in an axis 
system fixed to the vortex. In this system, the vortex 
is fixed and the airfoil, to which the grid is attached, 
moves past it. In this system 


Gi = 0 

Ifr = — Vr > 

= —( J 'v , 


(3.25) 


therefore, 

AG = A/ a Gr (3.26) 
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mile r ' 01 !. 3 A eliminates the need to com- 
pute G explicitly and hence the need to track a 

branch cut,. Equations (.3.20) through (3.26) are 

used to implement the split-potential method in the 
present algorithm. 

4 . Computational Results 

T J) e presentation of computed results is divided 
into four sections. Section 4.2 reports results ob- 
tained for simple flows (with no vortex interaction) 

Ihh°» V°, eStal ? l l Sh the Validity of the basic algo- 
rithm. Section 4.3 is an introduction to the results 

obtained in a HSI computation. The major features 
of the interaction are presented for a typical case 
the graphical presentation of data which describes 
the interaction is presented and explained. The sen- 
sitivity of the algorithm to both time step and ran- 
dom disturbances is discussed. Section 4.4 presents a 
comparison of the results obtained with the four vor- 
tex models. Comparisons with experimental data are 
a f S °, " la< Action 4.5 presents a parametric studv 
of the effects of some key parameters such as Mach 
number vortex strength, miss distance, vortex core 
size, and angle of attack. 

4.1. Executing the Algorithm 

4-1.1. Executing Baste Airfoil Solution 



1Z. 


..... . xunui,i re * sl(l,JHl convergence history for NA( - A 

0012 airfoil at several Mach mnnhers and o — ()°. 


A basic airfoil solution consists of the stead v flow 
around an a.rfo.1 at, a fixed angle of attack." The 
-steady-state version of the algorithm (eq. (2.59)) is 
used to generate this solution. Figure 12 depicts 

NAcTmir C ° nVerge,,Ce histor >' ^ an 

number Tt ' 3 ° = °° several Mach 

numbers. The residual is reduced by oscillating the 

I^eudo time At between the values 0.001 and 5 for 
the first 20 iteration steps. This oscillation drops the 
ratue ot R IIlax approximately 1 order of magnitude 
Ihe pseudo time is then fixed at At = 0.001 to 
minimize shock oscillation problems for the higher 
Mach numbers. Fixing At results in the unufual 
flattening of the curve at about 100 time steps After 

0001 iS ,,0 ‘ thp 0 P‘i*™m choice 
to reduce R. However, the current scheme provides 

suitable method for reducing R to an acceptable 
value over the widest, range of flow conditions. 

4-1 -2. Executing HSI Problem 

The solution of a blade- vortex interaction prob- 
lem proceeds in two steps. The first step is the com- 
putation of a .steady-state solution. The steady-state 


J J,,uu,ueu wun °r Without, the pres- 
ence of the vortex: however, for convenience the vor- 
tex is included. The vortex is initially fixed at some 
distance upstream (usually about 10 chords) and the 
steady flow about the airfoil is computed. After the 
s early-state solution is obtained, the vortex is al 
owed to move along a path determined bv equa- 
tions (3.4a). A specified path is produced by usine 
equations (3.4b) to compute vortex velocities "As the 
vortex moves through the computational grid the 
( .me step is varied by using the following relationship 
" uch was established by numerical experimentation: 


A( = f OAcospUy-f- 10.405)1 


10.405 


+ 0 



(4.1) 

Equation (4.1) is designed to provide a sufficient 
number of steps within each grid cell to assure proper 
resolution of the vortex. Thus, when the vortex is far 
from the airfoil At is large to minimize the number of 

a rS| r T n'r SOl ! ,,iQn: and when * near the 
loil, At is small to enhance accuracy. Application 

ot a varying time step reduces the time requirements 
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for computation by as much as a factor of 6 over a 
constant time step. 

4J.3. Data Recovery 

The velocities in the Bow field ate computed by 
the cent ral difference equations (2.33) and (2.34). 
the physical frame they become 

t,= + (4.2) 

v = + &vVy J 

Pressure is calculated by using the isentropic flow 
relations which lead to 

= (pT _ i) (4.3) 

' (l/2)pUi 

Lift on the airfoil is calculated by using the trape- 
zoid rule to integrate the pressures to get 


= 2 /’ C Ac,, 

J 0 


where c = sfxl The integration with respect to 
the square rod of the surface coordinateuK^ 
the accuracy in the leading-edge region. The Mach 

number is ^ 2 

ir + 1 


Ml = 


The density is computed directly from the Bernoulli 
equation (eq. (2.2)) 


(> = P* 1 + 


. - 2*1 - 


#5 - *1 


4.2. Results From Basic Airfoil Solution 

Before presenting the results for the HSI problem, 
it is useful to demonstrate the validity of the basic 
algorithm by comparisons with the results obtmned 
bv other researchers and by experiment for steady 
and unsteady flow problems. 

4 2 . 1 . Comparison With Original Algorithm 

The basic scheme used in the present method 
is an extension of a method developed by Steger 
and Caradonna (ref. 2). The first step m verifying 
the present method is to demonstrate the abi l 
to reproduce their results. The pre dieted mrdclmrd 
values of c p of a parabolic arc airfoil as it thicken. 
then thins for a free-stream Mach number of 0.8o 


Current method 

□ Reference 2 


-- 1 o 5 10 15 20 25 30 

Distance traveled, chords 

Figure 13. Pressure coellicients at nudd.ord for parabolic 
arc airlbil as it thickens then thins lor - 0.8o and 

r = 0.104. Small disturbance grid used. 

and maximum thickness r of 0.104 are P^ed m 
figure 13 The results for the present method were 
predicted by using a “small-disturbance type mesh 
(see step 1 in section 2.2); reference 2 also uses such a 
mesh Figure 13 demonstrates very close agreement 
between the present method and reference 2. 

Submit cal Flow 

Although HSI is, in general, a transonic flow prob- 
lem the demonstrated ability to compute subcritical 
flows is necessary to establish confidence in the al- 
gorithm. Figure 14 presents computed results for an 
NACA 0012 airfoil at, M x = 0.63 and a - 2° along 
with the results of Holst (ref. 25). 

Figure 15 presents calculated results 
by Michel Costes of the Office National d Etudes et 
de Recherches Aerospatiale (ONERA), w 10 la een 
provided with the present code) for Aejo»pat,ale 
RA16SC1 airfoil at M*. — 0.30 and o ■ 

suits are compared with experimental data generated 
h v ONERA and provided for this comparison, 
comparisons demonstrate the ability o ie 
algorithm to predict accurately the subcritical Ho 
about, an airfoil. The correlation with t he Aerospa- 
tiale RA16SC1 airfoil is particularly good considering 
the large difference between the surface grid line (ob- 
tained from the Joukowski airfoil) used to represent 
the airfoil surface and the actual airfoil surface. 


4.2.3. Supercritical Case 

A demonstrated ability to compute supercritical 
flow is necessary in order to proceed with the HSI 
problem. Figure 16 presents computed results tor 
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figure 15. Pressure coefficients for Aerospatiale RA16SC1 

airfoil at M x = 0.30 and a = 0°. Data from Michel 
C'ostes, ONERA. 


an NACA 0012 airfoil at A/ x = 0.8 and « = 0°. 
The surface pressures are compared with results 
generated by Bridgeman, Steger, and Caradonna 
(ref. 10), which were provided by him for this com- 
parison. Close examination of the figure reveals a 


slight difference in shock location between the two 
models. This difference is probably caused by the dif- 
ference in grid metric computation. The Bridgeman 
method uses the free-stream subtraction method to 
remove the metric truncation error and the present 
method uses consistent metric differencing. 

Figure 17 presents calculated results for an Aero- 
spatiale RA16SC1 airfoil at. A/ x = ().7f> aiK ] 0 — yo 
(provided by Costes). The results are compared 
with experimental data generated by ONERA. The 
correlation is good, except for a slight difference 
in shock location which is probably due to viscous 
effects on the real airfoil. 

4-2.4- Oscillating Airfoil 

The ability to predict unsteady flows is necessary 
in order to proceed with the HS1 problem. Figure 18 
presents calculated results for an NACA 0012 air- 
foil at My. = 0.755 and o = 0° and experimental 
results from Goorjian and Guruswamy (ref. 2(5). 
Comparison erf the two results shows excellent agree- 
ment, Figure 19 presents the results for an oscil- 
lating NACA 0012 airfoil. Figure 18 is the steady 
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Figure 16. Pressure coelfkienls lor NA( A 0012 airloil at 
= 0.80 hikI n = 0°. 

Exp. Calc. Sur. 


□ Upper 

O Lower 



Figure 17. Pressure coefficients for Aerospatiale RA1GSC1 
airloil at A/ x = 0.76 and <> = 0°. 

starting point . The airfoil is oscillat ing at a reduced 
frequency of k = ufc/l’x — 0.1(528 about o — 0 
with Omax = ±2.51. The predicted pressures are 
slightly different from the measured values because 
of the difficulty in exactly matching the test points. 
These results indicate that the algorithm can predict 
unsteady transonic flow. 

Figure 20 presents calculated results for an 
Aerospatiale RAK5SC1 airfoil with an oscillating 
25-percent flap at. = 0.30 and a = 0°. The re- 
sults are compared with experimental data for three 



Figure 18. Pressure coefficients for NACA 0012 airloil at 
= 0.755 and o = 0°. 

periods of motion. The |>ressures are presented in a 
format based on the following equat ion: 

= (16) 
k 

where (r,,)o is the steady pressure and <p is the phase 
angle of the oscillation. The correlation is good 
except for the flap pressures themselves. However, 
the grid is sparse in this region and the difference 
could easily be due to poor resolution. 

4.3. Introduction to HSI 

JJ.l. Generic Vortex Interaction 

Before presenting the vortex model comparison, 
a discussion of the basic features of the blade-vortex 
interaction will be useful. Figure 21 is a representa- 
tive plot, of lift coefficient versus vortex location for 
an NACA 0012 airfoil at a = 0° and = 0.60, and 
a vortex with an equivalent lift of q ,, = 0.400, and a 
vertical miss distance of 0.251 chord moving along a 
fixed path at constant speed. The computation uses 
the split-potential method to model the vortex. The 
several points of interest which occur during the in- 
teraction are labeled in the figure. Point 1 is the ini- 
tial condition of the airfoil in relatively undistorted 
flow (the vortex is far upstream). The sign of the 
vortex is such that, it, generates a downwash on the 
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Figure 19. Pressure coefficients for oscillating NACA 0012 airfoil at A/ x = 0.755 


airfoil, which has the same effect as reducing the an- 
gle of attack. When the vortex approaches, the pres- 
sure on the lower side of the airfoil diminishes and 
drives the lift coefficient to increasingly negative val- 
ues. This region is labeled 2 in the figure and is called 
the approach phase of the interaction. As the vor- 
tex approaches the leading edge of the airfoil, the lift 
rapidly diminishes to a minimum point, labeled 3. As 
the vortex passes under the airfoil the velocities in- 
duced by it become an upwash on a larger and larger 
portion of the airfoil. While the vortex is passing 
under the airfoil, the pressures change rapidly and 
the lift rebounds to positive values, labeled 4 ; this is 


the interaction phase. As the vortex approaches the 
trailing edge, the lift, begins to level out, labeled 5 
in the curve; this is called the recovery phase. Af- 
ter the vortex passes the trailing edge of the airfoil, 
the lift levels out and remains relatively constant as 
the vortex moves away. This region is labeled 6 in 
the figure and is called the departure phase. An in- 
teresting phenomenon occurs during the departure 
phase. The lift is found to rebound to a level which 
is higher than the initial condition and to return to 
i} 1 , 6 condition only after a considerable time. 

This behavior has been seen in most (if not all) other 
B\ I computational results reported in the literature 
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Figure 20. Pressure coefficients lor Aerospatiale RA1GSC1 
airfoil with oscillating 25-percent flap. A/x = 0.30; 
o=0°. 

but has not been explained yet. This decay effect 
may be explained by recalling that the HS1 problem 
is a limiting case of an infinite line vortex interact- 
ing with an infinite aspect ratio wing. As the vortex 
passes under the wing, it generates waves all along 
it. These signals cannot arrive simultaneously at any 
point. Therefore a “delay" effect is caused by the 
time lag in the arrival of waves from the rest of the 
wing. This “lag" could account for the delay in lilt 
recovery. This response characteristic is familiar m 


the field of acoustics but is unusual in aerodynamic 
problems. 

Figure 22 is a composite of the airfoil lift history 
and surface pressures at selected points during the 
interaction. The airfoil lift history is presented in 
the upper left-hand plot and the surface pressure 
coefficients are presented at points 1 through 5 in the 
other plots. These plots show that the presence of the 
vortex affects the airfoil primarily through changes 
in the lower surface pressure. The passage of the 
vortex is seen in the distortion of the lower surface 
pressure near the vortex location. (See plots for 2, 

3, and 4.) The plot, for 4 show’s a slight breakdown 
in the Kutta condition when the vortex is just past 
the airfoil trailing edge; this is shown later to be an 
effect of the time step. 

Figure 23 shows the time history of the lower 
surface pressures during BVI. Pressure coefficients 
as a function of time are presented at. 10 locations 
on the airfoil. The j>oint.s 2, 3, and 4, placed here 
for reference, correspond to the pressure plots in 
figure 22. Ftom this figure, it can be seen that the 
primary effect, is in the leading-edge pressure with a 
secondary effect propagating with the location of the 
vortex. The format, used in figures 21, 22. and 23 is 
used in the rest, of this report. 

Effects of Time Terms 

Because of the time linearization, testing the al- 
gorithm for sensitivity to the size of the time step is 
necessary. Recall that the time step is varied with 
vortex location (eq. (4.1)); this variation tends to 
increase the error when the vortex is far from the 
airfoil. Figure 24(a) shows the effect on airfoil lift 
history of reducing At by a factor of 2 from the val- 
ues computed bv equation (4.1). The figure shows 
the interaction of an NACA 0012 airfoil with a vor- 
tex of strength c/ ,. = 0.496, a miss distance of —0.433 
chord, and a free-stream Mach number of 0,536 ( i.e., 
the same condition as the sample case), figure 24(b) 
shows the interaction of an NACA 0012 airfoil with 
a vortex of strength c l v = 0.40, a miss distance of 
—0.251, and a free-stream Mach number of 0.80. The 
split- potential method was used to make these com- 
putations. Further reduction has no effect on the 
solution. The primary effect of reduced At is to 
improve the pressure coefficient prediction when the 
vortex is near the trailing edge (recall plot for 3 in 
fig. 22); that, is, the lower time step maintains a more 
accurate Kutta condition. This effect is a relatively 
minor one. Equation (4.1) was used for the following 
comparisons. Since the purpose' of these runs was 
to compare vortex models and explore parametric 
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Figure 24. Variation of lower surface pressure coefficient with 
vortex location lor typical BY I. 


effects, use of t he larger time . step makes little 
difference. 

However the reduced time step com put at ions do 
indicate the presence of “waves” superimposed on 
the basic solution especially for the departure phase 
of the high-speed case, dhese waves are generated 
as the vortex moves past nodes in the grid causing 
an abrupt change in local velocity. The wave ap- 
pears to behave exactly as an ordinary acoustic wave 
and is propagated with a Doppler effect. That is, 
a low-speed disturbance travels with equal strength 
upstream and downstream; a high-speed distur- 
bance travels with increased strength upstream. Fig- 
ure 25(a) is a plot of pressure at a fixed point, in 
space as the vortex passes at low speed; figure 25(b), 
at high speed. The high-speed disturbance is seen 
to occui after the vortex has passed the low-speed 
disturbance, before and after the vortex has passed. 
Reducing the time step in figure 24 more accurately 
captures this effect. These waves are seen to be mi- 
nor disturbances on the basic solution and do not 
affect the comparisons. 



(a) SuluTiUcal flow; a = 0.15; a = o°; Jl/ - 0 53G- 
y, — —0.443; q , — 0.406; NAOA 0012 airfoil; fixed path. 



(I>) Supercritical flow; a = 0.05; o = 0°; Jl/ = 0.80* 
U* ~ —0.251; C{ , ^ 0.40; NACA 0012 airfoil; fixed path. 


Figure 24. Effect of time step on airfoil fill during BYT 


One of the key unanswered questions associated 
with the use of the spit-potential model is the effect 
on the solution of temporal difference terms in G, 
( ^ iL,AG). (See eq. (3.14).) These terms can have 
a serious impact on a 3-D computation because of the 
geometric complexity of the wake and the difficulty 
of computing G ( for each wake element. Figure 26(a) 
shows the effects of the G/ terms on the integrated 
lift curve for a subcritical flow condition. From 
the figure, little difference is seen between the two 
curves. However, this flow condition is for low speed 
and the vortex is not very close to the airfoil (a so- 
called weak interaction). Figure 26(b) shows a sim- 
ilar comparison for a close supercritical interaction 
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G ( terms 




(l>) A/x = 0.80. 

Figure 25. Effect of'error wave propagation during BVI. 


(a “strong* interaction). Here the effect of the terms 
is much more apparent. In particular, the G( terms 
show a marked effect during the interaction phase 
and the departure phase. Notice especially the re- 
gion in which the slope of the curve undergoes a rapid 
change. This region corresponds to the vortex pass- 
ing through the airfoil shock. The combination of 
t hese velocities (shock and vortex) causes high pres- 
sures on the airfoil surface, which are reflected in the 

loading curve. 

In summary, the G, terms seem to have little 
effect on the solution except for the strong interaction 
cases. For the parametric studies, these terms are 
dropped frail any further computations involving the 
split- potential method in order to minimize run time 

and costs. 



(a) Suhcrit icaJ flow; a = 0.15; o - 0°; A/ x - 0.536; 
y _ -0.433: c/ , = 0.496; NACA 0012 airfoil; fixed path. 



(b) Supercritical flow'; a — 0.05; a — 0 , Ai x 0.80, 
y t - _o.25; cj , = 0.40; NACA 0012 airfoil; fixed path. 

Figure 26. Effect of G, terms on split-potential model. 

Jf.3.3. Comparison With Other Codes 

Figure 27(a) shows the airfoil lift history from 
the split-potential results and related computations 
of several other researchers (refs. 13, 15, and 17). 
The various other methods shown range from small 
disturbance algorithms to Euler methods. The in- 
teraction depicted is for an NACA 0012 airfoil at 
= 0.30 and a = 0° and for an equivalent lift 
of the vortex c iv = 0.40. The present method com- 
pares well with the results from the Euler equations 
(ref. 15), which has a more complete physical mod- 
eling, especially in prediction of minimum lift. Both 
the small disturbance (ref. 13) and the strongly im- 
plicit scheme (ref. 17) predict similar results. 
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(I>) Supercritical flow. A/ x =0.80. 

Figure 27. AirloU lilt variation witli vortex location Cor various algorithms. 


Figure 27(h) shows a similar comparison for 
the NACA 0012 airfoil at AI = 0.80 and a = 0°, 
c l,v = 0-40, and a vertical miss distance y,, of 
-0.25. These results were taken from a survey paper 


by Srinivasan and McCroskev (ref. 27). Figure 27 
demonstrates that the present method produces re- 
sults for the integrated lift, history comparable with 
those generated by other researchers. 


31 



Figure 28. Vortex-induced velocity on airfoil surlace. 


4.4. Comparison of Vortex Models 

4 . 4 .I. Vortex Velocity Field 

A straightforward way of comparing the four vor- 
tex models is to compare the induced velocity pro- 
duced on the airfoil by each model. Figure 28 
presents a series of vort ices located at the airfoil lead- 
ing edge and at selected vertical distances. These 
computations were made with a small disturbance 
mesh. (See section 2.2, step 1.) The computations 
were made for a flow field with no airfoil present in 
order to eliminate air foil -induced velocities from the 
results. The lifting-surface model is not indicated 
in the figure because it has the same effect, as the 
split-potential method. The angle-of-attack method 
in sharp contrast to the other methods is a simple 
step change in velocity. The branch-cut and sp it- 
potential methods both predict an impulsive type 


velocity curve. The sharpness of the impulse is mod- 
ified in the split-potential method by the vortex core. 
The velocity induced by the branch cut. is generated 
implicitly by the model. These velocities tend to be 
higher than the others because of the interaction of 
the four subvortices with the dense grid region near 
the airfoil. Distance for the airfoil reduces the size 
of the veloc ities and reduces the differences produced 
by the vortex on the airfoil. 

4.4.2. Related Experiment 

Comparison of the four vortex modeling meth- 
ods is made with the help of experimental data. 
Caradonna, Laub, and Tung (ref. 28) presented ex- 
perimental results for a rotor interacting with a vor- 
tex. The rotor had two blades with a constant 0012 
airfoil section. The blades were untwisted and the 
rotor had a teetering hub. The rotor aspect ratio 
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was /. A vortex was generated upstream of the ro- 
tor by a fixed, constant-section, NACA 0015 wing. 
Figure 29 depots the experimental setup. When the 
rotor blade is at an azimuth angle of 180°, HSI oc- 
curs. Pressure at, 10 locations on the airfoil surface 
were measured and these data were presented as a 

chord) 011 ° f time (V ° rteX IOCaU ° n) and s P ace (airfoil 

The data were collected for several rotor tip 
speeds and vortex locations. Two of these condi- 
tions are used as reference data in comparing the 
vortex modek The first condition is for a subcritical 
ow, - 0.536, and a vertical miss distance of 
w T1,e ^cond case is for a critical flow, 

J x ' ° and the same vertical miss distance. 

There are two factors that compromise the cor- 
re at i on. The first factor is the low aspect ratio of 
the rotor. Because the present method is strictly 
two-dimensional, one should expect to see a higher 
pressure predicted due simply to aspect ratio ef- 
fects. The second factor is the rotational velocity 
° the rotor Tl,e measurements used in the com- 
parison were made at an azimuth location of 180°. 
The effect, of the vortex on the blade, however be- 
gins much earlier. The blade is, therefore experi- 
encing a steadily changing free-stream flow modified 
by he vortex. The variable free stream is not mod- 
eed by the current method. For the subcritical flow 
this does not. pose a problem because unsteady ef- 
fects are small at low Mach numbers. However, the 
critical flow is much more sensitive to this effect. 

ie rotor is experiencing a high transonic speed at 
the azimuth location of 90° which decreases as the 
b ade moves forward. As the speed decreases, the 
shocks on the airfoil surface begin to collapse The 
vortex is encountered during this collapsing process 
Hie computation of this flow field requires a three- 
dimensional model complete with an accurate un- 
steady shock model; this is beyond the capability of 
the current method. J 

4-4-3- Sub critical Interaction 

Comparison of the various vortex models is made 
firs for a subcritical flow. The condition selected 
is the same as that used in the sample computation 
presented m the introduction, that is, an NACA 0012 
airfoil at o = 0°, M x = 0.536. a vortex strength 
/, U.4J6, and a constant miss distance of —0.433 

c ior .7 his condition should be relatively insensitive 
to unsteady effects and free from shock waves. 

Figure .10 shows lift versus vortex location for 
the four vort ex models. The four methods show lit - 
t e difference in the initial condition solution. As 


NACA 00 15 wing 



Section A-A 

Figure 29. Experimental measurement oil, I ade- vortex 
interaction. 


Method 



Figure 30. Variation of airfoil lil t with vortex location for the 
our vortex models. Subcritical flow; a = 0 . 15 ; a = 0°* 

A/*, zz 0.536; y, = -0.433; cy.,. = 0.496; NACA 0012 
airfoil; fixed path. 
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the vortex approaches, the curves begin to sepa- 
rate, reflecting the effect of the various methods 
on the flow-field solution. When the vortex is be- 
tween 2 and 3 chords upstream, the separation of 
the curves becomes large. This point can there- 
fore be considered the outer boundary of the close 
interaction. The curves continue to separate until 
they reach a maximum difference when the vortex is 
at the airfoil leading edge. During the interaction 
phase, the methods display several interesting fea- 
tures. The angle-of-attac k method predicts a rapi 
(almost instantaneous) change in lift as the vortex 
passes the quarter-chord because the airfoil under- 
goes an abrupt change in angle of attack from neg- 
ative to positive at this point; in effect the curve is 
inverted. An interesting feature of the branch-cut 
method is that during the interaction phase it pre- 
dicts a small “spike” in the loading curve. This spike 
occurs when the vortex is ill the densest part of the 
grid and reflects the locally high velocities predicted 
bv this method. This spike increases in size as the 
Mach number increases. For a very strong i» er ^' 
tion (A/ = 0.80, » P = -0.25, and c l v = 0.40). this 

spike leads to an instability in the algorithm which 
destrovs the solution, The split-potential and lifting- 
surface methods are in agreement overall even t hough 
they are somewhat different during the closest part o 
the interaction. Neit her method shows any unusual 
features in their predicted loading histories. 

Figure 31* presents pressure time histories for 
the four methods. Here the difference between the 
lifting-surface and split-potential methods is more 
apparent. The spike in the branch-cut loading curve 
is also more apparent in figure 31(c) than in the 
integrated data. Figure 31(a) highlights the change 
in pressure on the airfoil caused by the sharp change 
in lift predicted by the angle -of- at tack method. 

Figure 32 presents measured and calculated data 
for each of the methods. Comparisons can only be 
considered qualitatively valid because of the aspect 
ratio and unsteady rotational flow field of the ex- 
periment. The split-potential, lifting-surface and 
branch-cut methods all show qualitat ively good com- 
parisons but differ in detail. (Compare plots for 1 
and 2 in figs. 32(b), (c), and (d).) The angle-of- 
attack method is clearly not accurate lor this con- 
dition. (See plots for 1 and 2 in fig. 32(a).) This 
inaccuracy is caused by the sharp change m lift pre- 
dicted by the model which clearly does not occur in 
the experiment . 


•Figures 31 through 47 are at the end of section 4. 


4 . 4 . 4 . Cnttcal Interaction 


In section 4.4.3, the various vortex models were 
compared for a subcritical flow condition. A more 
interesting comparison is for a flow condition jus >e- 
low critical, that is, a flow condition which if undis- 
turbed would remain subcritical. The introduction 
of a vortex in such a flow field would be expected to 
drive the flow into a supercritical state. An NAt A 
0012 airfoil at a = 0° and M = 0.714 experiences 
this type of flow. 

Figure 33 presents lift versus vortex location for 
each of the models at this “critical” flow condition. 
Inspection of this figure shows qualitatively the same 
results as the subcritical case: the models begin to 
separate between 2 and 3 chords upstream, the angle- 
of-attack method predicts a sharp ‘lift inversion as 
before, the branch-cut method predicts a spike in u 
loading curve of increased size at this higher Mach 
number (0.714), and the split-potential and luting- 
surface methods both predict smooth curves. 


Figure 34 presents the pressure time histones on 
the airfoil. The angle of attack and lifting-surface 
methods both predict pressure histones which are 
similar to those at the lower Mach number. How- 
ever the branch-cut and split-potential curves are 
markedly different. Both methods indicate the pres- 
ence of shock waves although t he branch-c ut method 
is obsc ured by the presence of the spike. I he pres- 
ence of the shock wave in t hese models is more easily 
seen when comparing the measured and calculated 
data (fig. 35). The split-potential method predicts a 
shock which is also present 111 the data. (See plots 
for 2, 3, and 4 in fig. 35(d).) The branch-cut method 
predicts a sharp pressure peak which appears to be 
the beginning of a shock formation; however, the 
velocity spike may be delaying the formation The 
lifting-surface and angle-of-attack methods do not 
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The surface specification models also predict that 
the vortex affects the pressure on both surfaces of the 
airfoil. This is caused by the assumption of infinite 
signal speed which is inherent in these methods, lhe 
explicit methods show an effect on the lower surface 
only due to the additional time required lor the signal 
to arrive at the upper surface. 


^4.5. Summation 

The previous discussions show that all four mod- 
els produce qualitatively similar results for the inte- 
grated load but differ considerably in detail, lhe 
split-potential and lifting-surface methods produce 
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the most accurate curves for integrated loads. Di- 
rect. comparison between computation and experi- 
ment is not possible except in the most qualitative 
sense because the experiment contains large three- 
dimensional and unsteady effects which are not mod- 
eled in the theory. The angle-of-attack and branch- 
cut methods both generate spurious spikes in the 
loading curve. Comparison of the pressure time his- 
tories of the models shows considerable differences in 
the details of the pressure loading. Comparison of 
the methods for a critical flow condition shows that 
only the explicit models predict the presence of su- 
percritical flow. Only the explicit models are capa- 
ble of predicting the variation in local Mach number 
which is necessary to capture accurately the nonlin- 
ear nature of the flow field. Numerous computational 
experiments have indicated that the split-potential 
method is the most robust of these two methods. It. 
also allows for more control of t he vortex modeling 
because it. models the vortex as a. specified velocity 
field. In fact the method is not restricted to modeling 
vortices and can be used to predict the effect, of any 
flow field which is described by a known irrotationaj 
Potential function. Because of this versatility and 
robustness, the split-potential method is the recom- 
mended method. The rest of the calculated results 
presented in this report were generated by using the 
split-potential method. 

4.5. Parametric Sweeps 

This section contains a study of the effect of 
several key parameters on HSI. The split, -potential 
method is used in all comparisons. 

4-5.1. Effect of Mach Number 

One of the key parameters is Mach number. For 
an airfoil in steady flow, the effect, of increased Mach 
number is the increase in local velocit ies on the sur- 
face which leads to lower surface pressures. The vari- 
ation of the local pressure and lift coefficient can be 
predicted with good accuracy by using the Prandtl- 
Glauert correction. As the Mach number increases 
the local velocities increase until, finally, a shock 
develops on the surface of the airfoil. The linear 
I randtl-GIauert correction breaks down completely 
with the onset of local sonic flow, and the problem 
becomes nonlinear. 

The presence of a vortex in the flow simply adds 
an extra component to the local velocity because of 
the induced velocity field which the vortex generates. 
This extia velocity tends to induce the supercritical 
flow at a free-stream Mach number lower than the 
critical point for the airfoil alone. Figure 30 shows 
the effect that Mach number has on the blade-vortex 


interaction. The three curves are for an NACA 0012 
airfoil at o = 0° and q ,, = 0.40. The miss distance 
is -0.251 chord, the vortex core size is 0.05 chord, 
and the path of the vortex is specified bv setting 
the vortex velocity equal to the free-stream value. 
The Mach numbers are 0.4, 0.6, and 0.8. Almost no 
difference is seen between the subcrit.ical flows except 
during the interaction and departure phases. The 
supercritical flow, however, shows a broader pulse 
width than the subcriticaJ flows. This effect was also 
noted experimentally in reference 28 and is caused 
by the effect of signal propagation speed. The signal 
speed is the sum of the speed of sound and the flow 
velocity. The waves propagating from an upstream 
vortex arrive at the airfoil sooner for a higher speed 
interaction. A wave generated by a downstream 
vortex takes longer to arrive at the airfoil, and hence, 
the effect, of the vortex decreases more gradually if 
the propagation speed is increased. This is reflected 
in the figure as a separation of t he curves during the 
interaction and departure phases. Figure 37 presents 
t he pressure distributions on the airfoil surface for 
selected vortex locations. The curves reflect, the effect 
of Mach number and strongly resemble the variation 
one might, see for an isolated airfoil at negat ive angle 
of attack except, for the suction peaks induced when 
the vortex is very near to the blade surface. 

4-5.2. Effect of Vortex Strength 

A helicopter rota- airfoil interacts with a series 
of vortices generated by the preceding blade. The 
strength of these vortices is dependent upon the con- 
ditions under which the generating blade sheds them. 
Modeling vortex strength and understanding its ef- 
fect. are, therefore, important. The primary effect of 
vortex strength is to increase the peak velocities gen- 
erated by the vortex and consequently the velocities 
imposed on the airfoil. This effect can be deduced 
from the tangential velocity equation (eq. (3.3)). An 
interaction of particular interest would be the “criti- 
cal” interaction described in section 4.5.1. Figure 38 
is the integrated lift, history of an NACA 0012 air- 
foil for three different, vortex strengths (r/ ,, = 0.20, 
0.40, and 0.60) for the critical Mach number 0.714. 
The loading curve shows what appears to be a lin- 
ear variation with vortex strength. However the air- 
foil pressures (fig. 39) show an increasingly stronger 
shock wave on the airfoil surface. (See plots for 3 
and 4 in fig. 39.) The shock location and strength 
vary linearly with the vortex strength. For strong 
vortices (c t ,, = 0.60), the shock extends into the flow 
field and interacts with the vortex directly to cause 
very high suction peaks. (See plot for 3 in fig. 39(c).) 
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4.5.3. Effect of Miss Distance 

The location of t he vortices also affect the airfoil 
response. Consider again equation (3.3) and note 
that the geometry effect is primarily in the denom- 
inator. The effect of miss distance should be the 
inverse of the effect of q ,, . Figure 40 shows the ef- 
fect of miss distance on the integrated lift history of 
an NACA 0012 airfoil at M x = 0.714, o = 0°, and 
q ,, = 0.40. The miss distance varies from y v = -0.75 
to -0.25. The miss-distance variation in figure 40 is 
strictly a vertical displacement. Since the peak vor- 
tex velocity depends on radial distance, the effect of 
a vertical displacement does not become large until 
the proportion of r which it contributes is large. This 
effect begins to occur at approximately 1 chord up- 
stream from the airfoil leading edge. After this point, 
the change in loading appears to vary linearly with 
the inverse of the miss distance. The pressure dis- 
tribution curves (fig. 41) show a similar trend with 
respect to shock strengt h and location. 

4.5.4- Effect of Core Size 

The effect of core size on the airfoil response is 
more intricate. In equation (3.3), the vortex core 
size has the effect of modifying the denominator. As 
a increases, the denominator grows for a given value 
of ?*. Furthermore, the minimum value of lift occurs 
at the minimum value of the denominator when the 
vortex is at the leading edge. After this point, r 
increases and V decreases. At the point where r is a 
minimum, the value of a determines 1 . The effect 
of a is, therefore, to change the effect ive location of 
the vortex, and this is what gives rise to the phase 
shift, in figure 42. Figure 42 shows the pressure on 
the airfoil surface at selected vortex locations. As 
before, the shock strength and location vary directly 
with the inverse of a. 


4.5.5. Effect of Angle of Attack 


A helicopter rotor blade experiences a cyclic 
change in angle of attack as it rotates. The effect 
of the airfoil angje of attack on BVI is therefore im- 
portant- The angle of attack of the airfoil affects the 
response primarily by changing the initial conditions 
of the solution. Figure 44 shows the effect of a on 
the interaction between an NACA 0012 airfoil and 
a vortex of strength q t , = 0.40. The Mach number 
is 0.60 and a = 0°, 0.5°, and 1°. The lift history 
for the curve for a — 0° has been subtracted from 
the others to illustrate the difference between the in- 
teractions. As can be seen in the figure, the effect 
of angle of attack is purely linear for this suberiti- 
cal interaction. The airfoil pressure distributions are 
presented in figure 45. One interesting feature is that 
the suction peaks on the airfoil are reduced as the an- 
gle of attack increases. (See plots for 2 and 3.) I he 
positive angle of attack tends to offset, the effect, of 
the vortex which induces a negative lift. 

The effect of angle of attack at a. transonic Mach 
number is similar to that at subsonic conditions. 
Figure 46 presents integrated lift curves for a Mach 
number of 0.8 and the same conditions as figure 44. 
The pressure coefficients indicate a decrease in shock 
strength as the airfoil angle of attack increases. An 
inverse effect, would lead to increased shock strength 
with increasing negative angle of att ack. A combina- 
tion of high negative vortex strength and high neg- 
ative angle of attack would produce a very strong 
shock which could possibly lead to a shock-induced 
separation. The very high suction peaks (q, = -1.5) 
seen in the airfoil pressure distributions (plots for 3 
in figs. 47(a), (b), and (c)) occur as the vortex moves 
through the supersonic region near the airfoil surface. 
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(c) Branch-cut method. 
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(d) Split-potential method. 

Figure 31. Variation of airfoil pressure with vortex location. Snl, critical How. 
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Hgiire 34. Variation of airfoil pressure with vortex location at critical (low. 
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Figure 42. Filed of core size on variation of airfoil lift, o - 0° ; Moc - 0./ 14; y r 
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5* Concluding Remarks 

A study of the full-potential modeling of a blade- 
vortex interaction has been made. A primary goal 
of this study was to investigate the effectiveness of 
the various methods of modeling the vortex. The 
problem was studied within the context of a two- 
dimensional model problem, which represents one of 
the limiting types of blade- vortex interactions. The 
model problem restricts the interaction to that of an 
infinite wing with an infinite line vortex moving par- 
allel to its leading edge. This problem provides a 
convenient testing ground for the various methods 
of modeling the vortex while retaining the essential 
physics of the full three-dimensional (3-D) interac- 
tion. The flow field is assumed to l>e inviscid, irrota- 
tional, unsteady and, in general, transonic. 

A full- potential algorithm specifically tailored to 
solve the blade-vortex interaction (BY I) was devel- 
oped to solve this problem. The algorithm is based 
on a model first presented by Steger and Caradoiina 
(NASA TM-8121 1, AVRADCOM TR80-A-14). The 
algorit hm makes use of the unsteady mass conserva- 
tion and Bernoulli equations to form a full -potential 
model of the flow field. The system of equations is 
reduced to one equation by using a Taylor-senes ex- 
pansion of the temporal derivative of the density term 
in the conservation equation. The spatial derivatives 
are recast in delta form, with the density written at 
the previous t ime step. The stability of the algorit hm 
in transonic flow is assured through the use of upwind 
biasing oft he density in the flux terms. The flux met- 
rics are computed by the consistent metric method, 
which has been found to be superior to the so-called 
free-stream subtraction method that has difficulties 
with grid singularities. The equation is approxi- 
mately factored into convenient geometric parts in 
order to reduce the matrix to a compact form. A 
tridiagonal matrix inversion is used to solve for the 
updated potent ial solution. The model has the capa- 
bility to predict the steady and unsteady flow about 
an airfoil under subcritical and transonic flow condi- 
tions. ( bm pari sons of the results predicted are made 
witli those presented by other researchers and with 
experimental data. The comparisons indicate that 
the algorit hm is able to predict basic unsteady tran- 
sonic flow about an airfoil. 

The basic algorithm has been modified to include 
the effect of a vortex passing near the airfoil, four 
different methods of modeling of the vortex were 
used : 

1. The angle -of- at tack method 

2. The lifting-surface method 


3. The branch -cut method 

4. The split-potential method 

The angle-of-attack method uses the velocity field of 
a point vortex to compute a vortex -induced velocity 
at the airfoil quarter-chord. This velocity is then 
used to compute an effective angle of attack of the 
airfoil. This method is identical to techniques that 
are currently in use in comprehensive helicopter rotor 
analyses. The lifting- surface method is an extension 
of the angle-of-attack method in which the vortex- 
induced velocity is a function of chord wise distance 
on the airfoil surface. The branch-cut method is 
a flow-field vortex representation that makes use 
of a surface of potential discontinuity, the edge of 
which constitutes the vortex location. The effect 
of the vortex is implemented by imposing special 
differencing methods along the cut. In the split- 
potential method, the velocity field is split between a 
known field (induced by the vortex) and an unknown 
perturbation field due the airfoil. 

A side-by-side comparison of the four models was 
conducted. These comparisons included: 

1. Comparing generated velocity fields 

2. A subcritical interaction 

3. A critical interaction 

The subcritical and critical interactions are also com- 
pared with experimentally generated results. 

The following conclusions have been reached as a 
result of these comparisons: 

1. All the methods give qualitatively similar re- 
sults for integrated loads but differ consider- 
ably in the details. 

2. An artificial core is necessary in order to re- 
move the singularity in the vortex velocity 
equation. 

3. Only the explicit models predict the presence 
of shock waves for the critical interaction. 

4. The branch-cut method shows a strong sensi- 
tivity to the mesh configuration which leads 
to spurious waves in the solution, especially 
for transonic flow conditions. 

5. The angle-of-attack method is sensitive to miss 
distance and predicts a spike in the loading 
curve which is not present in the other meth- 
ods. This spike is caused by the abrupt, change 
in angle of attack as the vortex passes the air- 
foil quarter-chord. 

6. The lifting-surface method compares well with 
the split-potential method in computing in- 
tegrated loads, especially for subcritical flow. 
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However, the details of the pressure time histo- 
nes during the closest part of the interaction 
difTer significantly. Nevertheless this method 
has proven to be useful in 3-D problems to 
specify far -field vortex effects. 

7. The split-potential method is the most versa- 
tile and robust method. 

8. The temporal gradient terms in the split- 
potential method must be retained in order 
to predict strong interactions. This retention 
may cause a slight increase in computing t ime 
for 3-D problems, but this increase can be off- 
set. by using the lifting-surface method for vor- 
tices far from the blade. 

The split-potential model was used to make a sur- 
vey of some of the more crit ical parameters which 
affect the BVI. I he survey studies general flow pa- 
rameters such as free-stream Mach number and air- 
foil angle of attack and vortex parameters such as 
strength, core size, and miss distance. The results 
were computed at subcritical and supercritical free- 
stream Mach numbers. For the vortex parameters, 
the free-stream Mach number was chosen to be just 
subcritical to study the effect of the vortex on the 
formation of critical flow on the airfoil. Based on 
the survey results, the following conclusions were 
readied: 

1. Free-stream Mach number has little direct ef- 
fect on the basic BVI except to modify the sig- 
nal propagation speed which tends to broaden 
the loading signature. This effect was also 
noted experimentally in NASA TM-86005, 
I'SAAYSCOM TM-84-A4). However, the vor- 
tex can induce shocks on an airfoil which is in 
near-critical flow. 

2. The BVI loading response varies linearly 
with the vortex strength. The vortex induces 
shocks on the airfoil surface and the shock 
strength varies linearly with the vortex 
strength. For strong vortices, the shock may 
interact with the vortex directly to cause high 
peak pressures on the airfoil. 

3. The BVI loading response varies linearly with 
the inverse of the vortex miss distance. Vortex- 
induced shock strength also varies linearly 
with the inverse of the miss distance. 


4. The BVI loading response' varies directly with 
the inverse of vortex core size a. Furthermore, 
the core size shifts the point of minimum lift in 
direct proportion to a and the shock strength 
varies directly with a. 

5. The effect of airfoil angle of attack a is to 
shift the magnitude of the loading curve by 
a constant value which varies linearly with cv. 

A primary goal of this report has been to study 
the effect of the vortex model on the computation 
of the BVI problem. The selection of an appro- 
priate technique of modeling the vortex would be 
based on this study; this has been accomplished. 
The split-potential model has proven to be' the most 
versatile and robust, method of the currently avail- 
able techniques. The lifting-surface method has been 
shown to be a useful approximation to the split- 
potential method especially for far-field vortex spec- 
ifications. The next, logical step in this study is 
to extend these results to the 3-D rotor problem. 
In the process of accomplishing this, the method 
should be coupled with an existing comprehensive 
helicopter method, similar to that used by St r awn 
and Caradoiuia (AIA A-8fi-0079) . The 3-D model 
should include a. completed vortex wake model that 
uses a combined split-potential and lifting-surface 
method. Another interesting application of the split- 
potential technique would lie to use it to model lin- 
ear portions of the flow field (e.g., the rotational flow 
field). 

The modeling of the BVI continues to be a key 
problem in helicopter aerodynamics because it is a 
major determinant of vibratory loading and noise. 
This report has analyzed the interaction in two di- 
mensions and used this model problem to find the 
best means of determining BVI loading within the 
context of a finite-difference computation. The ex- 
tension to three dimensions should build directly on 
this work. Beyond this point, the greatest, problem 
will be to find an efficient and accurat e way to predict, 
the three-dimensional struct ure of the rotor wake. 


NASA Langley Research ('enter 
Hampton, VA 23681-0001 
April 22, 1997 
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Appendix A 

Computational Grid 

The computational grid used in the present 
method is based on the streamlines and potential 
lines around a Joukowski airfoil in incompressible 
flow and at. a = 0°. The intersection of this set of 
lines forms an orthogonal H-type mesh. This typeo 
mesh is usefiil for three reasons. 

1 ](. closely conforms to the shape of an airfoil, 

thereby accuracy is increased 

2. It is orthogonal which simplifies coding 

3. The gridlines align with the free stream away 
from the airfoil 

The solution of the Joukowski airfoil is a classical 
problem in aerodynamics. The solution is based on 
the conformal mapping technique which arises from 
complex variable theory. The details of this well- 
known derivation are provided for students of the art. 

The stream function '!'(£.»/) and the potential 
function $(£,»f) are related through the differential 
equations 

<)& rW 
Ox ~ Oy 

0$ _ (W 
dy Ox 

Both functions are solutions to Laplace s equation 

V 2 # = V 2 '!) = 0 (A'2) 

Therefore, the stream function and the potential 
function may be combined into an analytic function 
of a complex variable r*: 

«;(;*)=$ + »♦ (A3) 

The funct ion u’(- + ) is often referred to as the complex 
potential. When given a known function «’(-*)■ the 
potential field may be deduced by setting 4> equal to 
the real part of w(:*) ■ 

The solution of Laplace’s equation for problems 
with complicated boundaries is achieved with the aid 
of conformal mapping. The problem is mapped from 
the c* plane to another plane in which the bound- 
aries are simplified. The primary restriction on the 
transformation is that the mapping function be an- 
alytic. Joukowski presented a function which trans- 
forms the flow about a circle to that about an airfoil 
shape. (See ref. 29.) From this solution, obtaining 



an equation for the streamlines and potential lines 
about this airfoil shape is possible. 

Consider the flow field about a circle of radius b, 
then the complex potential is 


w = f'x I z + — 


where 


- = £ + ir l 


From equation (A3), 


$ + w = Utk. £ + >y + 


£ + iy 


Equating real and imaginary parts gives 


* = 1 + TTrTi 


6- 

£ 2 + >i‘ 


* = 1 “ TTTdl 


lr 

£ 2 + y 2 


The velocities and "F ,, are 


= U 


(£ 2 + n 2 ) 2 + b'Hy 2 - £ 2 ) 


(£ 2 + y 2 ) 2 


.. b 2 i^e>i 
(£ 2 + y 2 ) 2 

The derivatives ^ r, and are 


= 1 


, (£ 2 + y 2 ) 2 - b 2 {y 2 - £ 2 ) 


(£ 2 + y 2 f 


, .. >> 2 V^iy 

1 (£ 2 + y 2 f J 

The slope of the streamlines is 

dy^Z±L = (A10) 

d£ *£ (£ 2 + y 1 ) 2 - b 2 (£ 2 - y 2 ) 

Equation (A 10) may be integrated to produce 

?/;+! = C*A£,- + yt (All) 

where (* is equal to Equation (All) is for a line 
of constant stream funct ion. 


74 



The grid is generated with the following steps: 

1 Produce a satisfactory stretched Cartesian 
grid by any method 

2. Use the (£, tj) coordinates along the front of 
the mesh as a start to i ntegrate equation ( A ]()) 
until you reach the aft face of the mesh (this 
solves for the streamlines around the circle) 

3. Transform the circle solution by using the 
Joukowski transformation to produce the air- 
foil solution 


4. Select an appropriate distribution of point 
along the airfoil streamline 

5. Interpolate to find the potential at each of 
these points 

6. Find the location on each of the off-airfoil 
streamlines which have matching potential 
values 

7. Form the new mesh with the resulting set of 
points 
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Appendix B 

Derivations of Boundary Conditions 


The equation for Of is 



<;* 



(B4) 


There are two special equations which are used 
in implementing the boundary conditions in the full- 
potential algorithm. Both these equations make use 
of special applications of the Bernoulli equation. 

The first equation is used along the aft face of 
the computational grid to ensure that p — 1. A 
special equation is required because of the presence 
of branch cuts which specify a jump in potential 
and hence preclude the use of a Dirichlet condition. 
The equation is used to determine the velocity which 
is used in the flux computation. For p = 1, the 
Bernoulli equation yields 

2#r + U$f +V#,, = M'ic (Bl) 


Equation (B4) is used in the flux computation 
along the aft. face of the comput ational mesh in place 
of the usual backward difference of *. 

The second equation is used to compute the con- 
vection of T downstream of the airfoil trailing edge. 
Since the wake cut is a shear layer, density and V 
are constant across it. The jump in potential is 

r= ||*||= ( Br >) 

Using the Bernoulli equation gives 

2 ($ l)r + A| [ df (*/)]' 2 = 2 dr (*/+ n 

+ A,[^(#/ + r)] 2 (B6) 


The aft face of the computat ional grid is far from the 
airfoil (80 chords). It is therefore valid to assume 


V = 0 

At =(i 

‘he = o>- 4 4 / x c 


(B2) 


With these assumpt ions obtain 


Mr + JWxW( = ( B3 ) 


Eliminating terms gives 

2r r T 2A[*fl f T A) rj; = 0 ( B 7) 

The higher order term is dropped for convenience to 
give 

IV+(U)r { = 0 (B8) 

The term (U) is simply the average of the velocities 
above and below the cut. Equation (B8) is used to 
determine the converted value of 1 along the wake 
cut. 
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Appendix C 

Time Linearization Operator 

The conventional form of an operator which is 
derived from the Bernoulli equation and acts on the 

difference *' v+1 - $ A is given by equation (2.22), 
which is 




(JL v£\ a ’ 

\dr + + Or] ) 


The operator arises from the linearization of the den- 
sity which is necessary to maintain strong conserva- 
tive form. The linearization takes the form 


Op m 

p- Po + (<-u 

The subscript o represents a neighboring known state 
or solution. If ($ - $ 0 ) is small, for example 
(* A+1 - * A ) or (#,- - #,_]), where i = nAt and 
* = the error due to expanding p is second- 

order accurate and is no greater than that usually 


made in time differencing. Equation (Cl) is obtained 
by using the Bernoulli equation 


1 = p» + [l + yp ( M i - 2*r - U* f - V$„)j 


+ 


7 - 1 


— t y-rl-n* 


i-V+l 


- 


- 2 A,** (**+'- **)- 2 A 3 *^ 


;V + 1 - * V 
’) V 


(C2) 


Combining terms gives 


(> = (>» 



( o_ ud va \ 

V dr + d£ + (h) ) 


so that the operator becomes 


(#- v+1 - v ) 
(C 3 ) 



vW- + — +— Y 

\<>t ^ ot (h) ) 
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